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The  draw  resonance  instability  of  a bicomponent  coextrusion  fiber  spinning  has 
been  studied  theoretically  and  experimentally.  For  the  theoretical  analysis  a simple  model 
has  been  adopted  in  which  a Newtonian  and  an  upper-conveeted  Maxwell  fluid  is  the 
core  and  the  skin  layer,  respectively.  This  model  is  chosen  to  investigate  a two-phase 
flow  in  which  the  two  fluids  have  quite  different  extensional  rheology.  The  results 
indicate  that  the  viscoelastic  skin  layer  has  a stabilizing  influence  in  which  the  onset  of 
draw  resonance  of  a Newtonian  fluid  is  delayed  to  a higher  draw  ratio  by  the  coextrusion. 
An  unexpected  behavior  of  the  coextrusion  flow  is  predicted  in  that  the  highest  draw 
ratio  for  a stable  coextrusion  spinning  can  be  larger  than  that  of  the  single-phase  spinning 
of  either  fluid  (Newtonian  or  upper-convected  Maxwell)  under  certain  conditions.  In  the 
experimental  study  a linear  low-density  polyethylene  and  a low-density  polyethylene 
were  used  as  the  core  and  the  skin,  respectively.  These  materials  have  been  chosen  to 
simulate  the  flow  situation  considered  in  the  theoretical  investigation.  While 


xi 


the  theoretical  model  is  too  simple  to  provide  a quantitative  agreement  with  the 
experiment,  various  aspects  of  the  experimental  observations  appear  to  be  in  qualitative 
agreement  with  the  theoretical  predictions. 

A novel  processing  method  has  been  developed  for  the  fabrication  of  plastic 
optical  fibers  in  which  a modified  coextrusion  fiber  spinning  technique  is  used.  Lab  scale 
experiments  have  demonstrated  that  the  new  method  is  capable  of  fabricating  defect-free 
clad  fibers  using  several  different  types  of  materials.  The  new  method  is  possible  because 
the  mechanics  of  the  bicomponent  spinning  flow  is  controlled  by  the  viscoelastic 
cladding  material  under  typical  spinning  conditions  as  the  process  modeling  study  has 
indicated. 


CHAPTER  1 
INTRODUCTION 


1.1  Description  of  Fiber  Spinning 

An  important  industrial  process  for  manufacturing  artificial  fibers  is  called  the 
fiber  spinning  process  in  which  a polymer  material  is  continuously  extruded  downward 
through  a spinneret  or  through  a die  to  form  an  extrudate  of  circular  (or  noncircular) 
cross  section.  Somewhere  downstream  of  the  spinneret  the  extrudate  is  pulled  by  a take- 
up  machine  to  form  a thin  filament.  The  technique  of  spinning  synthetic  fibers  has  made 
significant  progress  during  the  past  thirty  years.  New  spinning  processes,  as  well  as  new 
varieties  of  man-made  fibers,  have  been  developed.  For  a given  fiber-forming  material, 
different  spinning  techniques  can  sometimes  produce  fibers  having  markedly  improved 
physical  and/or  mechanical  properties.  The  fiber  industry  has  therefore  made  continuous 
efforts  to  modify  existing  processes  and  develop  new  ones.  There  have  been  two  recent 
developments  in  fiber  spinning  technology:  the  production  of  "shaped"  fibers  and  the 
production  of  "coextrusion"  (or  bicomponent)  fibers.  A shaped  fiber  is  one  whose  cross 
section  is  noncircular.  A coextrusion  fiber  consists  of  two  different  components,  extruded 
side  by  side.  Each  of  them  has  its  advantages  over  fibers  made  in  the  conventional  way. 

One  of  the  conventional  types  of  spinning  processes  is  commonly  known  as  "melt 
spinning"  in  which  a fiber  can  be  made  by  extruding  molten  polymer  and  allowing  the 
thread  to  solidify.  In  melt  spinning,  fluid  elements  are  subjected  to  a strong  elongational 
flow.  Unlike  shearing  flows,  the  elongational  flows  generally  constitute  a ffee- 
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boundary  problem.  The  elongational  flow  of  a fiber  spinning  process  has  been  given 
substantial  attention  in  the  last  three  decades,  due  to  the  increasing  importance  of  the 
economics  of  synthetic  fiber  production  and  the  quality  of  the  product.  The  typical 
characteristic  of  the  elongational  flows  including  the  melt  spinning  can  be  found  in  the 
monographs  by  Ziabicki  ( 1 976),  Petrie  ( 1 979),  and  Denn  ( 1 980). 

The  melt  spinning  process  is  shown  schematically  in  Figure  1.  The  fiber-forming 
zone  consists  of  three  regions:  die  swell  region,  draw-down  region  and  solid  region.  Die 
swell  is  caused  by  the  stress  relaxation  as  the  material  leaves  the  confining  geometry  of 
die  or  spinneret.  Viscoelasticity  of  polymers  typically  results  in  a much  larger  die  swell 
than  purely  viscous  materials.  The  dependence  of  the  die  swell  on  rheological  properties 
and  on  the  die  geometry  is  quite  complicated,  and  the  theoretical  prediction  of  the  die 
swell  is  as  yet  incomplete.  Many  experimental  results  have  shown  that  for  a given 
polymer  at  a given  temperature  die  swelling  increases  as  output  rate  increases,  and  as  die 
length  and  die  diameter  decrease  (Bergonzoni  and  DiCresce,  1966).  The  die  swell  region 
is  followed  by  the  draw-down  region  in  which  the  molten  polymer  is  drawn  to  a fine 
filament  as  the  axial  velocity  increases  to  the  predetermined  take-up  velocity.  Because  of 
the  die  swell  the  fiber  diameter  near  the  spinneret  can  not  be  predicted  accurately. 
However,  this  region  of  uncertainty  is  restricted  to  a small  region  within  a few  spinneret 
diameters  of  the  die  exit,  and  most  of  the  spinning  length  corresponds  to  the  draw-down 
region.  In  the  analyses  of  the  fiber  spinning  flow,  therefore,  the  initial  filament  diameter 
is  typically  taken  as  either  the  maximum  observed  value  or  the  spinneret  diameter  itself. 
The  location  of  the  point  of  solidification  may  be  controlled  precisely  by  rapid  quenching 
with  blowing  air,  a water  bath  or  a chill  roll. 

At  the  early  stage  of  the  analyses,  researchers  usually  have  neglected  any 
interactions  between  the  fiber  and  the  surrounding  medium  (i.e.,  air).  Under  these 
assumptions  the  flow  is  isothermal  and  no  shear  or  normal  stresses  act  on  the  filament 
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Figure  1 . Schematic  of  fiber  spinning  process 
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boundary  (Matovich  and  Pearson,  1969;  Shah  and  Pearson,  1972a,  1972b;  Spearot  and 
Metzner,  1972;  Donnelly  and  Weinberger,  1975;  Denn,  Petrie  and  Avenas,  1975;  Fisher 
and  Denn,  1975a,  1975b,  1976;  Chang  and  Denn,  1979,  1980;  Chang,  Denn  and  Geyling, 
1981;  Kase  and  Nakajima,  1980;  Larson,  1983).  A nonisothermal  analysis  may  simulate 
the  real  situation  more  closely.  However,  the  inclusion  of  energy  equation  coupled  with 
the  equation  of  motion  only  makes  the  system  to  be  more  complicated  without  seriously 
affecting  the  fundamental  characteristics  of  the  fiber  spinning  flow  (Kase  and  Matsuo, 
1965,  1967;  Shah  and  Pearson,  1972a,  1972b;  Nakamura,  et  al.,  1972,  1974;  Fisher  and 
Denn,  1977;  Gupta  and  Metzner,  1982). 

The  most  important  aspect  of  the  fiber  spinning  flow  is  that  the  fiber  is  slender 
and  its  diameter  varying  slowly  in  the  flow  direction.  This  allows  an  asymptotic  analysis 
which  result  in  a one-dimensional  representation  of  the  fiber  spinning  flow  (Schultz  and 
Davis,  1982,  1984a;  Park,  1990). 

1.2  Definition  of  Draw  Resonance 

When  polymer  melts  are  subjected  to  fast  drawing  processes  such  as  fiber 
spinning,  extrusion  coating,  or  film  casting,  an  instability  phenomenon  called  "draw 
resonance"  is  often  encountered  in  which  the  fiber  radius  (or  the  film  thickness)  and  the 
take-up  force  fluctuate  periodically.  This  instability  occurs  when  the  draw  ratio,  which  is 
defined  as  the  ratio  of  the  take-up  to  the  extrusion  velocity,  is  greater  than  a certain 
critical  value  and  the  fluctuation  is  observed  to  be  sustained  with  a well-defined  period 
and  amplitude.  The  terminology  of  "draw  resonance"  was  coined  by  Christensen  (pp.752, 
1962):  "Surging,  or  draw  resonance,  is  defined  as  a non-uniformity  of  gauge,  or  coating 
weight,  in  the  machine  direction."  and  its  periodic  nature  was  first  described  by  Miller 
(PP- 137,  1963):  "Sustained  periodic  fluctuation  of  fiber  radius  (or  film  thickness)  which 
occurs  above  a certain  critical  take-up  velocity  (or  draw  ratio)."  The  draw  resonance  has 
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been  studied  extensively  ever  since  due  to  its  industrial  importance.  Among  the  various 
drawing  processes,  the  fiber  spinning  has  been  studied  most  extensively  probably 
because  its  geometry  is  simpler  than  the  other  processes.  Despite  their  geometric 
differences,  however,  the  development  and  understanding  of  the  fiber  spinning  process  is 
directly  applicable  to  extrusion  coating  and  film  casting  processes  since  the  mathematical 
descriptions  in  modeling  these  processes  are  analogous. 

The  theory  of  draw  resonance  was  first  developed  independently  by  Kase  and 
Matsuo  (1965,  1967)  and  by  Pearson  and  Matovich  (1969).  These  analyses  were  carried 
out  for  an  isothermal  spinning  of  a Newtonian  fluid  under  the  assumption  that  the  inertia, 
gravity,  and  surface  tension  effects  were  negligible.  Their  results  predicted  that  the 
critical  draw  ratio  for  the  onset  of  draw  resonance  is  20.21 . The  theory  was  then  extended 
to  power  law  fluids  and  viscoelastic  fluids  (Denn  et  al.,  1975;  Fisher  and  Denn,  1976; 
Phan-Thien,  1978,  1987;  Petrie,  1979;  Larson,  1983;  Phan-Thien  and  Caswell,  1986; 
Gupta  et  al.,  1986;  Mackay  and  Petrie,  1992).  These  studies  have  shown  that  the  shear 
thinning  behavior  of  polymeric  fluids  has  a destabilizing  effect  whereas  the 
viscoelasticity  has  a stabilizing  effect.  The  theoretical  predictions  are  generally  in 
reasonable  agreement  with  many  experimental  observations  at  least  qualitatively 
(Christensen,  1962;  Donnelly  and  Weinberger,  1975;  Matsumoto  and  Bogue,  1978; 
Chang  and  Denn,  1979;  Demay  and  Agassant,  1985).  A thorough  survey  of  numerous 
studies  on  the  subject  was  given  in  a review  paper  by  Petrie  and  Denn  (1976),  and  an 
updated  list  of  references  can  be  found  in  more  recent  papers  by  Denn  (1980,  1983)  and 
Mewis  and  Petrie  (1987). 


1.3  Instabilities  in  Fiber  Spinning  Flow 

Many  different  types  of  instabilities  are  encountered  in  polymer  processing.  From 
a practical  standpoint  it  is  important  to  identify  the  parameters  which  dictate  the 
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instability  and  to  understand  the  instability  behaviors  as  affected  by  the  parameters. 

In  the  fiber  spinning  process,  there  may  be  three  different  types  of  instabilities; 
melt  fracture,  draw  resonance  and  fiber  breakage.  The  melt  fracture  phenomenon,  which 
results  in  surface  roughness,  is  common  to  most  extrusion  processes.  This  instability 
occurs  if  an  excessive  shear  stress  is  imposed  on  the  polymer  flow.  Although  it  is 
generally  agreed  that  this  instability  may  be  related  to  the  stick-slip  phenomenon  at  the 
polymer-metal  interface  (i.e.,  die  surface),  theoretical  understanding  of  it  is  only  at  the 
seminal  stage  at  the  present  time.  Such  melt  fracture  often  limits  the  output  rate  of  a fiber 
spinning  system,  since  successful  spinning  usually  requires  a smooth  extrudate. 

The  other  types  of  spinning  instability  arise  not  from  the  shear  flow  in  the  die  but 
from  the  free-boundary  flow  in  the  drawing  region  between  the  spinneret  and  the  take-up 
point.  The  draw  resonance,  in  which  the  fiber  diameter  fluctuates  periodically,  may  be  a 
kinematic  instability  as  it  is  predicted  without  all  external  and  surface  forces  except  the 
viscous  force.  Only  the  polymers  with  a very  low  viscoelasticity  are  susceptible  to  this 
instability  since  a large  viscoelasticity  may  alter  the  flow  kinematics  significantly. 

The  polymers  with  a large  viscoelasticity  exhibit  fiber  breakage  instead  of  draw 
resonance.  Since  the  breakage  may  occur  as  the  stress  imposed  on  the  fiber  exceeds  the 
cohesive  strength  of  the  material,  it  may  not  be  an  instability  phenomenon  in  a fluid 
mechanical  sense. 

In  this  thesis,  the  major  emphasis  is  on  the  theoretical  and  experimental 
investigation  of  the  draw  resonance. 

1.4  Coextrusion 

All  the  analyses  described  in  section  1 .2  are  for  a "single-phase"  flow,  and  few 
studies  have  been  devoted  to  a "multi-phase"  flow.  A multi-phase  flow  refers  to  a flow 
situation  in  a coextrusion  process;  two  or  more  polymers  are  processed  simultaneously  to 


7 


produce  composite  materials  such  as  multilayer  films  or  cladded  fibers  (e.g.,  optical 
fibers).  The  coextrusion  process  has  been  known  for  many  years  and  has  been  broadly 
used  for  numerous  applications  in  the  plastics  industry.  Nevertheless,  there  exist 
relatively  few  fundamental  studies  on  the  subject  (Han,  1975;  Waters  and  Keeley,  1987). 
As  was  indicated  in  many  studies  on  single-phase  fiber  spinning,  the  characteristics  of 
draw  resonance  (i.e.,  critical  draw  ratio,  amplitude  and  period)  depend  on  the  fluid 
rheology.  Thus,  a unique  flow  situation  is  encountered  in  coextrusion  processes  when 
the  rheological  properties  of  the  constituent  polymers  are  quite  different  from  one 
another. 

One  such  example  may  be  the  coextrusion  of  a linear  low-density  polyethylene 
and  a low-density  polyethylene.  The  former  is  known  to  be  very  prone  to  draw 
resonance  whereas  the  latter  is  not.  The  spinnability  of  a low-density  polyethylene  is 
typically  limited  by  the  breakage  of  spinline  as  the  stress  in  the  melt  may  exceed  the 
cohesive  strength  of  the  material  (Ziabicki,  1976).  The  difference  in  the  stability 
behaviors  of  the  two  materials  may  result  from  the  fundamental  differences  in  their 
rheological  properties.  A low-density  polyethylene,  which  has  a high  degree  of  long 
chain  branching,  has  a larger  relaxation  time  than  a linear  low-density  polyethylene. 
Consequently,  the  former  tends  to  show  a distinct  strain  hardening  behavior  when 
subjected  to  an  elongational  flow  whereas  the  latter  does  not  (Chang  and  Lodge,  1972; 
Spearot  and  Metzner,  1972;  Minoshima  and  White,  1986a,  1986b).  If  these  two 
materials  are  coextruded,  their  flow  mechanics  may  be  strongly  influenced  by  the 
rheological  interactions  of  the  two  materials. 

1.5  Scope  of  this  Study 

This  thesis  is  a theoretical  and  experimental  study  of  the  fundamentals  of  a two- 
phase  coextrusion  fiber  spinning.  In  chapter  2 one-dimensional  equations  are  presented 
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for  the  steady  flow  of  a two-phase  coextrusion  fiber  spinning.  As  pointed  out  previously, 
the  one-dimensional  equations  are  obtained  as  a result  of  an  asymptotic  analysis  in  which 
the  aspect  ratio  of  the  slender  fiber  is  the  small  parameter.  Simple  constitutive  models 
(Newtonian  and  upper-convected  Maxwell  models)  are  chosen  to  represent  the  two-phase 
system  so  that  the  competing  influence  of  viscous  and  viscoelastic  forces  is  investigated. 
Some  discussions  and  computations  are  also  given  to  explain  the  boundary  conditions 
which  have  been  used  for  the  present  analysis. 

In  chapter  3 the  stability  of  the  steady  flow  to  an  infinitesimal  disturbance  is 
investigated.  The  stability  criteria  are  given  by  the  neutral-stability  curves  not  only  for 
the  single-phase  flow  (Newtonian  or  upper-convected  Maxwell  fluids)  but  also  for  the 
two-phase  flow  (the  combination  of  a Newtonian  and  a UCM  fluid).  Considering  the 
periodic  nature  of  the  draw  resonance,  the  linear  stability  analysis  is  useful  and 
informative. 

In  chapter  4 experimental  investigation  of  the  draw  resonance  is  presented.  A 
linear  low-density  polyethylene  and  a low-density  polyethylene  have  been  used  for  the 
study  as  the  core  and  the  skin  layer,  respectively.  This  choice  of  materials  has  been  made 
to  simulate  the  flow  situation  considered  in  the  theoretical  investigation.  The  theoretical 
model  is  too  simple  to  expect  a quantitative  agreement  with  the  experiment. 
Nevertheless,  various  aspects  of  the  experimental  observations  appear  to  be  in  qualitative 
agreement  with  the  predictions  of  the  linear  stability  analysis.  The  results  in  fact  show 
similarities  in  various  trends.  The  critical  draw  ratios,  the  ratios  of  fiber  maximum  to 
minimum  radius,  and  periods  are  measured  and  analyzed.  The  results  show  that  the  core- 
layer material  is  stabilized  when  the  skin-layer  material  is  a viscoelastic  material.  This 
observation  may  be  of  important  industrial  implication. 

In  chapter  5,  a modified  coextrusion  fiber  spinning  process  is  applied  to  fabricate 
plastic  optical  fibers  using  non-conventional  materials.  Initially  the  experiment  was 
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conducted  specifically  for  the  on-going  efforts  to  develop  radiation  hard  scintillating 
optical  fibers  for  the  Superconducting  Super  Collider.  This  new  approach  was  then 
extended  to  make  plastic  optical  fibers  using  various  polymeric  materials  which  cannot 
be  formed  into  optical  fiber  by  the  conventional  process  technique  of  thermal-drawing. 
As  the  major  emphasis  of  the  experiments  was  on  the  processing  aspects,  little  attention 
has  been  given  to  monitoring  the  optical  properties  of  the  fiber  thus  fabricated. 

The  summary  and  conclusion  of  this  dissertation  are  finally  given  in  chapter  6. 


CHAPTER  2 

STEADY  FLOW  ANALYSIS 
2.1  Introduction 

The  fiber  spinning  process  has  been  investigated  to  date  from  several  different 
aspects;  steady-state  characteristics,  stability  behaviors  and  the  die  swell  under  applied 
axial  force.  The  early  investigations  focused  on  the  single-phase  fiber  spinning  of 
isothermal  Newtonian  fluid  (Kase  and  Matsuo,  1965,  1967;  Matovich  and  Pearson, 
1969).  During  the  last  several  decades,  the  reaserchers  have  considered  the  power-law 
fluid  and  various  viscoelastic  models  (Shah  and  Pearson,  1972c;  Denn,  Petrie  and 
Avenas,  1975;  Ishihara  and  Kase.  1976;  Petrie,  1978).  The  essential  difficulty  in 
analyzing  the  motion  of  a viscoelastic  liquid  is  the  choice  of  the  constitutive  equation. 
The  constitutive  or  rheological  equation  of  state  is  the  relation  between  stress  and 
deformation  of  a material.  Petrie  (1979)  states:  "We  do  not  know  what  equations  we 
should  be  solving  but  are  interested  in  the  way  the  choice  of  constitutive  equation 
interacts  with  the  dynamics  and  affects  the  mathematical  problem  to  be  solved  (pp.2)." 

The  first  theoretical  analysis  on  the  two-phase  fiber  spinning  flow  was  conducted 
by  Park  (1990).  As  in  his  study,  the  present  analysis  also  adopts  a simple  system  in  which 
a Newtonian  and  an  upper-con vected  Maxwell  fluid  are  assumed  to  be  the  core  and  the 
skin  layer,  respectively.  This  simple  model  is  chosen  to  investigate  a rather  extreme 
situation  where  the  two  fluids  have  quite  different  rheological  properties.  It  is  assumed 
that  the  flow  is  isothermal  and  that  only  viscous  and  viscoelastic  forces  are  important 
thus  neglecting  the  influence  of  inertia,  gravity,  and  surface  tension. 
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The  present  analysis  is  limited  to  the  draw-down  region  where  the  fluids  are 
stretched  forming  a thin  fiber.  In  this  region  the  change  of  the  fiber  radius  with  the  flow 
direction  is  very  small.  Consequently,  one-dimensional  approximation  is  applied  at  the 
outset  in  which  the  axial  velocity,  pressure,  and  stress  fields  are  assumed  to  be 
independent  of  the  radial  coordinate.  The  one-dimensional  approximation  was  proven  to 
be  valid  in  the  draw-down  region  by  asymptotic  analyses  (Schultz  and  Davis,  1982;  Park, 
1990,  1991)  and  numerical  simulations  (Fisher,  Denn  and  Tanner,  1980;  Keunings, 
Crochet  and  Denn,  1983).  Throughout  the  analysis  the  shear  viscosity  ratio  of  the  core 
layer  and  skin  layer  is  assumed  to  be  unity  while  the  relaxation  time  of  the  skin  layer  as 
well  as  its  flow  rate  fraction  are  varied. 

An  asymptotic  analysis  for  the  steady  flow  under  the  prescribed  conditions  has 
been  given  for  the  case  in  which  the  applied  stress  (or  draw  force)  is  infinitely  large 
compared  to  the  elastic  shear  modulus  of  the  skin-layer  material  (Park,  1990).  The  result 
of  the  asymptotic  analysis  is  extended  here  to  include  the  steady  state  solutions  for  a 
finite  draw  force.  The  steady  state  solution  serves  as  a base  flow  for  the  linear  stability 
analysis  presented  in  the  next  chapter.  The  steady  state  solutions  indicate  that  the 
viscoelasticity  of  the  skin  layer  plays  a significant  role  in  determining  the  overall 
mechanics  of  the  flow.  Thus,  a strong  influence  of  the  viscoelasticity  may  be  also 
anticipated  on  the  stability  of  the  flow  even  when  the  flow  rate  of  the  skin  layer  may  be 
smaller  than  that  of  the  Newtonian  core  layer. 

2.2  Problem  Formulation  and  Scalings 

The  schematic  of  a bicomponent  fiber  spinning  flow  is  shown  in  Figure  2.  The 
core  layer  and  skin  layer  are  assumed  to  be  a Newtonian  and  an  upper-convected 
Maxwell  fluid,  respectively.  The  region  of  interest  for  this  analysis  is  limited  to  the  draw- 
down region  where  the  materials  are  stretched  under  the  applied  axial  tension.  The  one- 
dimensional approximation  is  adopted  at  the  outset  in  which  the  flow  is  assumed  to  be 
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Figure  2.  Schematic  of  coextrusion  fiber  spinning 
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isothermal  and  axisymmetric.  Under  this  assumption,  the  stress  and  axial  velocity  are 
uniform  over  the  cross-section  of  the  fiber  at  each  axial  locations.  This  approximation  is 
known  to  be  invalid  in  the  immediate  vicinity  of  the  die  exit  where  the  extrudate  swell 
occurs.  However,  the  extrudate  swell  is  restricted  to  a very  small  region  and  the  one- 
dimensional flow  is  quickly  approached  within  about  two  radii  from  die  exit.  Thus,  the 
extrudate  swell  region  may  be  excluded  from  the  analysis  without  causing  a serious 
problem.  Therefore,  the  origin  of  the  system  is  taken  to  be  approximately  at  the  point  just 
below  the  extrudate  swell  region  in  which  the  rate  of  change  of  fiber  diameter  with  the 
axial  distance  is  very  small  (dR/dz«l).  In  most  of  the  studies  on  melt  spinning  the  same 
approach  has  been  taken  regarding  the  origin  of  the  coordinate  system. 


2.2.1  Spinning  Equations 

Under  the  assumption  that  the  gravity,  inertia,  and  surface  tension  are  negligible 
compared  with  viscous  and  viscoelastic  forces,  the  continuity  equation  and  the 
momentum  equation  are  given  as  follows: 


dA 
— + 

dt 


d_ 

dz 


(Aw)  = 0 


+ (\-f)A{(xca  -O- 


(2.1) 

(2.2) 


Here  A is  the  cross-sectional  area  of  the  fiber,  w the  axial  velocity  in  the  flow  direction. 
The  T ..  represent  the  extra  stresses  in  the  radial  rr,  circumferential  66,  and  axial  zz 

directions,  and  the  superscripts  s and  c denote  the  skin  and  core  layer,  respectively.  The 
fraction  of  skin-layer  flow  rate  is  expressed  by  / The  two  terms  in  the  first  curly  bracket 
of  (2.2)  are  the  axial  tension  bom  by  the  upper-convected  Maxwell  skin  layer  while  the 
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two  terms  in  the  second  curly  bracket  of  (2.2)  are  the  axial  tension  in  the  Newtonian  core 
layer. 

2.2.2  Constitutive  Equations 

The  rheological  expressions  are  needed  to  describe  the  relationships  between 

stresses  and  deformation  for  the  skin  and  core  material.  For  a uniaxial  elongational  flow, 
the  velocity  gradient  tensor  (Vy)  and  the  rate-of-deformation  tensor  D are  expressed  as 


Here  the  D is  defined  by  X[(^T)  + (^T)r]-  Since  off-diagonal  elements  of  Vy  are  zero, 
Vy  = Z).  For  the  Newtonian  core  layer  and  the  upper-con vected  Maxwell  skin  layer,  the 

constitutive  equations  are  written  as 


1 dw 
2~dz 


0 0 

!*:  o 

2 dz 


(Vy)  = D=  0 


(2.3) 


0 


dz 


f = 2flcD 


(2.4) 


(2.5) 


where  \ is  the  relaxation  time  of  the  upper-convected  Maxwell  fluid;  |ic  and  |is  are  the 

shear  viscosities  of  the  Newtonian  and  the  upper-convected  Maxwell  fluid.  These 
material  parameters,  |ic,  and  }is,  are  constant  in  the  given  constitutive  models. 


15 


The  substitution  of  Eq.(2.3)  into  Eqs.(2.4)  and  (2.5)  gives  the  radial  rr, 
circumferential  88,  and  axial  zz  components  of  the  constitutive  equations  for  the  core 
layer  and  skin  layer  as 


, e dw  , dw 


dz 


dz 


XL  + A 


dx^rr  dx3  dw 


\ 


/ 


To#  A 


dr*  dtL  . dw 
— - + w— + 4 — 


V 


dt 


dz 


dz 


r 


t^+A 


V 


dx3  dx3^  dw 

—r ^ + W— —~2f77 

dt  dz  dz 


(2.6, 2.7, 2.8) 
(2.9) 

(2.10) 

(2.11) 


2.2.3  Scalings  and  Dimensionless  Governing  Equations 

The  scalings  for  the  nondimensionalization  are  as  follows: 


L F 

(r,z)~(R0,L),  t , h-~w0,  A~A0,  x 

wo  Ao 

where  R0  is  the  radius  of  the  fiber  at  z=0  and  L is  the  draw  span;  w0  and  A0  are  the 
average  axial  velocity  and  the  cross-sectional  area  at  the  origin  (z=0);  F is  the  axial  draw 
force  applied  at  the  take-up  (z=l).  Under  these  scalings,  the  dimensionless  continuity 
equation  remains  the  same  as  equation  (2.1)  if  the  variables  in  (2.1)  are  considered  to  be 
dimensionless.  Combining  the  Eqs  (2.2),  (2.6),  (2.7),  and  (2.8),  the  momemtum  equation 
is  given  in  a dimensionless  form  as 
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+ 3(\-f)GMA^ 

az 


(2.12) 


The  superscripts  s are  dropped  here  for  convience  as  all  stress  components  appearing  in 
the  equations  are  for  the  skin  layer. 

The  constitutive  equations  for  the  upper-convected  Maxwell  fluid  are  given  by 


r^+De 


dr 


71 


V 


dt 


/ 


rn  + De 


dr 


IT 


dt 


r m + De 


f dr 


68 


dt 


(2.13) 

(2.14) 

(2.15) 


Here  the  time  dependence  is  preserved  as  they  will  be  used  in  the  following  chapter. 

The  dimensionless  parameters  appearing  in  the  above  equations  are  defined  as 
follows: 

Deborah  number 

De  = ^- 
L 


Shear  viscosity  ratio 


Dimensionless  inverse  axial  force 


Skin-layer  flow  rate  fraction 


Here  Qs  and  QT  are  the  flow  rate  of  the  skin 


M = 


Ml 


Q _ M 4H0 
FL 


f = 


_ Qs 


Qt 


layer  and  the  total  flow  rate,  respectively. 
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A rigorous  derivation  of  the  one-dimensional  Equations  (2.1),  (2.12),  (2.13), 
(2.14)  and  (2.15)  can  be  found  in  Park  (1990)  in  which  the  method  of  asymptotic 
expansions  is  applied  in  the  small  parameter  e which  is  defined  as  the  ratio  of  the  radial 

to  axial  characteristic  length  scales  R0/L.  The  most  significant  assumption  in  deriving 
Equations  (2.1),  (2.12),  (2.13),  (2.14)  and  (2.15)  may  be  the  radial  uniformity  of  the 
extra  normal  stresses  in  the  skin  layer.  While  the  velocity  independence  of  the  radial 
coordinate  arises  naturally  from  the  asymptotic  analysis,  there  is  no  clear  indication  that 
the  extra  normal  stresses  in  the  skin  layer  are  radially  uniform.  If  the  initial  stress 
condition  at  z= 0 is  independent  of  r,  however,  the  stress  field  becomes  radially  uniform 
for  all  z.  Thus,  for  simplicity,  it  has  been  assumed  that  the  stress  field  of  the  skin  layer  is 
also  independent  of  the  radial  coordinate.  This  assumption  may  not  be  a severe  restriction 
considering  the  successes  of  many  theoretical  analyses  built  on  the  same  assumption  for 
single-phase  flows. 


2.3  Boundary  Conditions 

At  the  origin  of  the  draw-down  region,  the  cross-sectional  area  of  the  fiber  and 
the  axial  velocity  are  given  as  A0  and  H'0,  respectively.  The  take-up  velocity  and  the  draw 
force  cannot  both  be  specified.  If  the  take-up  velocity  is  fixed,  then  there  exists  a 
corresponding  force  which  is  necessary  to  obtain  the  prescribed  take-up  velocity.  Here 
the  applied  force  at  the  take-up  point  is  adjusted  so  that  the  desired  value  of  the  final 
fiber  radius  (or  the  axial  velocity  at  z=L)  is  obtained.  Thus,  the  boundary  conditions  of 
the  cross-sectional  area  and  the  velocity  of  the  fiber  are  described  as 


A* 

= 4 

at 

z = 0 

(2.16) 

w - 

-~wo 

at 

z = 0 

(2.17) 

w - 

ZWL 

at 

z = L 

(2.18) 
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After  nondimensionalization  using  the  scalings  described  in  section  2.2.3,  the  cross- 
sectional  area  and  the  velocity  conditions  at  the  both  ends  of  the  fiber  are  given  in 
dimensionless  form  as 


A = 1 

at 

z = 0 

(2.19) 

w=\ 

at 

z = 0 

(2.20) 

w = Dr 

at 

z = 1 

(2.21) 

Here  Dr  is  the  draw  ratio  defined  as  the  ratio  of  the  take-up  to  the  extrusion  velocity 
w(l)/  w(0)  (or  the  ratio  of  the  initial  to  the  final  cross-sectional  area  of  the  fiber). 

One  underlying  assumption  in  specifying  these  velocity  conditions  is  that  no-slip 
occurs  at  the  interface  between  the  two  fluids;  that  is,  the  velocity  is  uniform  over  the 
entire  cross-section  of  the  fiber.  The  normal  stresses,  on  the  other  hand,  are  discontinuous 
across  the  fluid-fluid  interface  although  they  are  radially  uniform  within  each  layer.  The 
two  terms  in  the  curly  bracket  of  Eq.(2. 1 2)  represent  the  axial  tension  bom  by  the  skin 
layer  whereas  the  last  term  in  (2.12)  is  the  axial  tension  in  the  Newtonian  core.  Unless 
the  draw  ratio  or  the  Deborah  number  is  very  small,  the  tension  in  the  viscoelastic  skin 
layer  is  typically  much  greater  than  that  in  the  Newtonian  core;  that  is,  most  of  the 
tension  on  the  spinline  is  supported  by  the  skin  layer. 

In  addition  to  the  velocity  conditions,  the  stress  state  at  z= 0 needs  to  be  specified 
for  the  viscoelastic  skin  layer.  The  initial  stress  condition  accounts  for  the  strain  history 
of  the  fluid  to  that  point,  and  we  assume  that  it  is  given  as 


r*  (0)  = 'o 

U°)='o-l 

Tee(°)=^(0) 


(2.22) 

(2.23) 

(2.24) 
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Here  the  value  of  /0  is  an  arbitrary  constant.  Although  the  numerical  assignment  for  this 
value  is  arbitrary,  numerical  calculations  have  shown  that  the  mechanics  of  the  flow  is 
rather  insensitive  to  the  value  of  t0  (see  section  2.6).  For  most  of  the  calculations 
presented  in  this  study,  numerical  value  of  10  was  assigned  to  t0. 

While  arbitrary  values  can  be  assigned  to  xn  and  xm  regardless  of  x„ , it  was 

chosen  that  xn  = t0  - 1 and  xm  = xn  for  all  the  results  throughout  this  thesis.  For  a single- 
phase flow  of  an  upper-convected  Maxwell  fluid  (i.e.,/=l),  specification  of  at  z=0 
happens  to  be  unnecessary  since  its  value  is  determined  to  be  T„  (0)  - 1 by  the  axial  force 
balance.  Thus,  the  choice  of  T„.(0)  = t0-1  for  the  two-phase  flow  corresponds  to  the 
situation  that  the  upper-convected  Maxwell  fluid  in  the  two-phase  flow  experiences  the 
same  strain  history  as  its  single-phase  flow.  It  may  be  a sound  choice  for  the  comparison 
between  single-  and  two-phase  flows.  As  pointed  out  previously,  however,  the  flow 
mechanics  is  not  seriously  affected  by  the  numerical  value  of  /0  or  by  any  other  choice  of 
T In  section  2.6,  the  negligble  influence  of  /0  is  discussed  in  detail. 

The  initial  value  of  xm  is  chosen  to  be  the  same  as  ^(O)  in  order  to  drop  the 
(t#,  - xn)  term  in  Equation  (2.12).  This  term  accounts  for  the  hoop  stress  resulting  from 

the  concentric  core-skin  structure  of  the  flow.  As  will  be  shown  in  section  2.5,  the 

spinning  flow  is  also  insensitive  to  the  initial  value  of  the  hoop  stress.  Thus,  the  given 

choice  again  does  not  seriously  affect  the  conclusions  of  the  present  study.  Furthermore, 
when  r„.(0)-  1^(0),  it  can  be  shown  from  (2.14)  and  (2.15)  that  xm  = x„  for  all  z.  Thus 

the  momentum  equation  (2. 12)  is  simplified  to 


d_ 

dz 


fA(ru-Xn)  + 3(\-f)GMA 


(2.25) 


This  equation  is  more  convenient  than  (2.12)  for  the  comparison  between  single- 
and  two-phase  flow.  If  f is  set  to  zero,  Equation  (2.25)  is  reduced  to  the  one  for  a 
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Newtonian  single-phase  flow.  If/=1,  on  the  other  hand,  the  well-known  equation  for  the 
spinning  of  an  upper-convected  Maxwell  fluid  is  recovered. 

The  system  of  Equations  (2.1),  (2. 13), (2. 14)  and  (2.25)  along  with  the  boundary 
conditions  (2.19)  through  (2.23)  is  to  be  solved  throughout  all  this  work.  The  specific 
choices  for  the  boundary  conditions  may  appear  to  restrict  the  flow  situation  to  a very 
specific  case.  Due  to  the  insensitivity  of  the  flow  mechanics  to  the  initial  stress  condition, 
however,  the  present  results  in  fact  represent  the  general  situation  of  the  two-phase  flow. 
Further  justification  of  the  initial  stress  condition  is  given  in  sections  2.5  and  2.6,  and  an 
unexpected  interesting  implication  of  (2.24)  is  also  discussed  in  section  2.5. 


2.4  Steady  State  Solutions 

At  steady  state  the  time  derivatives  d/dt  are  set  to  zero,  and  the  steady  system  to 
be  solved  becomes 


dz 


(AaW“)=0 


d_ 

dz 


fA*(TZ--C)  + X\-f)GMA 


ss 


dw“ 

dz 


= 0 


tZ+De 
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V 


dz 


dz 


J 


= 2 G 


dw 


ss 


dz 


C + De 


rr 


V 


dz 


dz 


J 


dz 


(2.26) 

(2.27) 


(2.28) 


(2.29) 


The  boundary  conditions  are 


vT(0)  = l 


(2.30) 
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H'"(°)  = l (2.31) 

w“(l)  = Dr  (2.32) 

O0)  = 'o  (2.33) 

C(0)  = ^o~1  (2.34) 


The  integration  of  (2.26)  and  (2.27)  results  in 


f5  -T“  = — 

" "■  / 

A simple  manipulation  of  (2.28)  and  (2.29)  using  (2.35)  and  (2.36)  yields  the  following 
set  of  nonlinear  differential  equations  for  numerical  integration: 


w“  -3(\- f)GM 


dw* 

dz 


(2.36) 


dw“ 

dz 


(2.37) 


du 


ss 
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ss 


dz  Z{\- f)GMDe  w 
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0 ~f)M 


f K 


\ 


3 G Dew “ Gw a 
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drZ 

zz  


dz  Dew 


ss 


[2m“(G+DO-t^] 


(2.38) 

(2.39) 


Here  superscripts  ss  denotes  the  steady  state,  and  the  new  variable  uss  represents  the 
derivative  dwss/dz. 

Using  the  boundary  conditions  at  z=0  (i.e.,  Equations  (2.31),  (2.33), (2.34)),  the 
initial  value  of  uss  (i.e.,  dwss/dz)  is  derived  from  (2.36)  as 
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dwu  _ 1 

dz  ~ 3GM 


at  z = 0 


(2.40) 


Using  the  boundary  conditions  (2.31),  (2.33)  and  (2.40),  Equations  (2.37),  (2.38)  and 
(2.39)  can  be  integrated  numerically  to  determine  w"  and  r“  for  specified  values  of/  M, 

De  and  t0.  Once  w"  and  are  determined,  A”  and  follow  from  (2.35)  and  (2.36). 
The  fourth  order  Runge-Kutta-Gill  method  was  used  for  the  integration. 

The  inverse  axial  force  G is  not  a prescribed  parameter.  Instead  the  draw  ratio  Dr 
is  specified  at  z=l  (Equation  (2.32)).  Thus,  the  value  of  G is  determined  by  a shooting 
method  so  that  the  boundary  condition  (2.32)  is  satisfied.  In  Figures  3,  the  draw  ratio  Dr 
is  given  as  a function  of  G for  various / (/N),  0.2,  0.5,  and  1.0)  as  M=\,  De= 0.02,  and 

fo~  1 0.  For  a Newtonian  single-phase  flow  (i.e.,/=0),  an  analytic  expression  is  known  for 
the  relationship  between  Dr  and  Gmv: 


Dr  = exp 


1 


\ 


\5G 


ox 


Gm 


1 

3 In  Dr 


(2.41) 


(2.42) 


Here  GNW  is  the  inverse  axial  force  for  the  Newtonian  fluid,  which  is  defined  as 
G = H / FX.  Using  the  definition  of  G for  the  two-phase  flow  given  in  page  16, 
Equation  (2.41)  can  be  rewrtiien  as 


Dr  = exp 


\ 

/ 


3 GM 


(2.43) 


23 


Thus,  when  the  stability  of  the  flow  is  not  considered,  we  may  say  that  a Newtonian  fluid 
can  be  stretched  indefinitely  as  the  draw  force  is  increased.  On  the  other  hand,  Upper- 
convected  Maxwell  fluid  cannot  be  stretched  indefinitely  since  the  draw  force  becomes 
infinity  beyond  a certain  draw  ratio.  The  maximum  obtainable  draw  ratio,  which  is  51  for 
De=  0.02,  is  predicted  analytically  by  the  asymptotic  theory  (see  Appendix  B).  In  a 
single-phase  flow  of  either  a Newtonian  or  an  upper-convected  Maxwell  (UCM)  fluid 
(/= 0.0  and  1.0,  respectively),  the  draw  ratio  increases  monotonically  with  increasing 
draw  force.  However,  the  relationship  between  Dr  and  G is  shown  to  be  not  necessarily 
monotonic  for  two-phase  flows  as  indicated  in  Figure  3.  For  the  prescribed  condition 
(i.e.,  M=  1,  £>e=0.02  and  1 0) , the  curve  for/O.5  follows  the  monotonic  trend  of  the 
single-phase  UCM  fluid  approaching  Dr=  49  as  G -»  0.  When  f=  0.2,  however,  Dr  reaches 
a maximum  value  of  about  71  (point  c of  Figure  3)  before  the  infinite-stress  asympotote 
(Dr=50)  is  approached.  Thus,  when  Dr  is  between  50  and  71,  there  exist  two  steady  state 
solutions  for/=0.2  (e.g.,  points  a and  b of  Figure  3 for  Dr=  60). 

The  existence  of  the  maximum  draw  ratio  at  a finite  G for  a small  value  of  / 
seems  physically  conceivable.  When  the  draw  ratio  is  small,  the  viscous  force  dominates 
the  flow  mechanics  whereas  the  viscoelastic  force  does  when  the  draw  ratio  is  large. 
Thus,  a transition  may  occur  from  a viscous  force  dominant  mode  to  a viscoelastic  force 
dominant  mode  as  the  draw  force  (or  the  draw  ratio)  is  increased.  For  a relatively  large 
value  of/  (e.g.,  /=0. 5 in  Figure  3),  the  transition  may  occur  at  a draw  ratio  which  is 
smaller  than  the  asymptotic  value  for  G— >0.  Consequently,  Dr  increases  monotonically 

with  decreasing  G.  If /is  small,  on  the  other  hand,  the  viscous  force  is  dominant  until  the 
draw  ratio  is  much  larger  than  the  value  predicted  by  the  asymptotic  theory  (/= 0.2  curve 
in  Figure  3).  As  the  draw  force  is  further  increased,  the  viscoelastic  force  becomes 
significant  acting  against  further  stretching  of  the  fluid  resulting  in  a maximum  Dr  at  a 
finite  value  of  G. 
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G 


Figure  3.  Draw  ratio  as  a function  of  inverse  axial  force  at  De= 0.02, 
M=  1,  /o=10;  a:  G=0.023,  Dr==  60;  b:  G=0.064,  Dr=  60;  c: 
G=0.046,D/=7 1 
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This  behavior  of  the  coextrusion  flow  has  an  interesting  implication;  that  is,  the 
highest  draw  ratio  for  a stable  spinning  can  be  larger  than  that  of  the  single-phase 
spinning  flow  of  either  fluid  if  the  region  be  in  Figure  3 happens  to  be  stable.  It  may  be 
noted  that  the  Newtonian  single-phase  flow  is  unstable  when  Z>>20.21  (Pearson  and 
Matovich,1969).  As  will  be  discussed  in  Chapter  3,  the  low-G  branch  of  the  multiple 
steady  solutions  (i.e.,  the  portion  of  the/=0.2  curve  where  G<0.05  including  point  a)  is 
shown  to  be  always  unstable.  The  high-G  branch  of  the  solution,  on  the  other  hand,  is 
either  stable  when  De  is  large  or  goes  through  two  transitions  with  increasing  Dr  when 
De  is  small.  Thus,  when  De  is  large  enough  so  that  the  high-G  branch  solution  is  stable 
without  the  two  transitions,  it  is  possible  to  obtain  a much  higher  draw  ratio  with 
coextrusion  than  obtainable  with  either  the  UCM  fluid  or  the  Newtonian  fluid. 

Similar  implication  may  be  expressed  using  a different  type  of  plot.  In  Figure  4a 
and  4b,  the  inverse  axial  force  G is  plotted  as  a function  of  Deborah  number  De  at  two 
different  draw  ratios  Dr= 20  and  30,  respectively.  The  scale  of  inverse  axial  force  G in 
Figure  4 is  normalized  by  the  value  for  the  Newtonian  single-phase  flow  case,  Eq.(2.42): 


G 


Gm 


3Gln  Dr 


(2.44) 


In  Figure  4,  a draw  ratio  of  20  is  chosen  for  an  illustration  purpose  because  it 
corresponds  to  the  limiting  draw  ratio  of  a Newtonian  fluid  for  stability  (Pearson  and 
Motovich,  1969;  Denn,  Petrie  and  Avenas,  1975). 

These  figures  represent  the  maximum  value  of  De  beyond  which  the  given  draw 
ratio  (Dr=  20  or  30)  is  unattainable.  In  Figure  4a  for  Dr=  20,  the  maximum  De  for _/=1.0  is 
0.046.  This  value  indicates  that  an  upper-convected  Maxwell  fluid  whose  De  is  larger 
than  0.046  can  not  be  drawn  beyond  Dr=  20.  Thus,  if  Dr  is  desired  to  be  20,  the  material 
time  constant  should  be  small  enough  so  that  De  is  smaller  than  0.046.  When  f=  0.5,  the 


3G(lnDr) 
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Figure  4.  Inverse  axial  force  as  a function  of  Deborah  number,  De : 
(a)  Dr= 20,  /0=  10. (asymptotic  solution,  G— > 0:  1.  De=0.046, 
2.  De=0.047,  3.  Z)e=0.050);  (b)  £>r= 30,  t0=  10. (asymptotic 
solution,  G— >0:  1.  De=0.031,  2.  Z)e=0.032,  3.  Z)e=0.033). 
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De 


Figure  4.  (continued). 
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maximum  De  for  Dr=  20  is  allowed  to  be  only  slightly  larger  (i.e.,  0.047).  When  f=  0.2, 
however,  De  can  be  as  high  as  0.061  which  is  substantially  larger  than  0.046.  Thus,  in  a 
coextrusion  flow,  an  upper-convected  Maxwell  fluid  can  be  drawn  to  a higher  draw  ratio 
than  with  a single-phase  flow.  Similar  behavior  is  described  in  Figure  4b  for  Dr=  30. 

An  analytic  expression  for  the  maximum  Deborah  number  at  a given  draw  ratio  is 
given  in  the  work  of  Park  ( 1 990): 
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(2.45) 


When/^1  and  f0=l,  this  expression  is  reduced  to  Demax=  1 /(Dr- 1 ) which  is  for  the  single- 
phase flow  of  UCM  fluid  as  derived  by  Denn,  Petrie  and  Avenas  (1975). 

The  computaional  results  of  the  present  analysis  are  shown  to  be  very  close  to  the 

prediction  of  the  analytic  expression  supporting  the  accuracy  of  the  present  numerical 
calculation. 

The  multi-valuedness  of  the  coextrusion  flow  at  J^0.2  represented  in  Figure  3 is 
also  manifested  in  Figures  4a  and  4b.  When  De  is  between  0.050  and  0.061,  two  steady 
solutions  exist  for  Dr=  20  and/=0.2;  that  is  the  lower  and  the  higher  G solutions  which  is 
similar  to  a and  b of  Figure  3.  In  Figure  5 a similar  plot  is  given  but  for  various  values  of 

Dr  when^is  fixed  at  0.2.  It  appears  that  the  multi-valuedness  persists  over  a broad  range 
of  Dr. 

In  Figure  6,  the  steady  state  velocity  profiles  are  given  as  a function  of  the  axial 
coordinate  z for  M=  1,  De=  0.02  and  /o=10  at  two  different  values  of  draw  ratio  (Dr=  30 


3G(lnDr) 
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Figure  5.  Inverse  axial  force  as  a function  of  Deborah  number,  De  at 
different  draw  ratios  when/=0.2,  /0=  1 0,  M=\. 


30 


Z 


Figure  6.  Steady  state  axial  velocity  profiles:  (a)  Dr=30,  De= 0.02, 
r0=  1 0,  M=  1;  (b)  Dr=  40,  De=  0.02,  /0=10,  A/=l. 
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Figure  6.  (continued). 
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and  40).  It  is  indicated  in  the  figures  that  the  velocity  profiles  for  large  f(f=  1.0  and  0.5) 
become  more  linear  with  increasing  Dr  as  the  viscoelasticity  effect  of  the  skin  layer 
become  more  significant.  According  to  the  theory  of  single-phase  flows,  the  UCM  fluid 
flow  is  stable  at  both  Dr=  30  and  40  whereas  the  Newtonian  flow  is  unstable.  Thus,  it 

may  be  of  interest  to  determine  whether  the  two-phase  flows  are  stable  or  not  at  these 
draw  ratios. 


The  difference  xm  - xn  in  Equation  (2.12)  accounts  for  the  hoop  stress  resulting 
from  the  concentric  core-skin  structure  of  the  two-phase  flow.  Although  xm  and  x 

satisfy  the  same  differential  equation  (2.14)  and  (2.15),  they  are  not  necessarily  the  same. 
However,  if  the  initial  conditions  for  x^  and  xn  are  the  same  or  if  the  isotropic  pressure 

in  the  skin  layer  is  assumed  to  be  radially  uniform,  the  hoop  stress  term  in  (2.12)  is 
dropped  since  either  condition  results  in  xm  = xn  for  all  z. 

In  the  presence  of  the  hoop  stress,  the  set  of  equations  to  be  solved  for  the  steady 
state  is  as  follows: 


2.5  Effect  of  Hoop  Stress 


az 


(2.26) 


(2.46) 


(2.28) 


(2.29) 
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(2.47) 


The  intergration  of  the  steady  equations  (2.26)  and  (2.46)  associated  with  the  boundary 
conditions  (2.30)-(2.34)  gives 


(2.35) 


w“-  3(1  -f)GM 


dw* 

dz 


(2.48) 


Using  (2.35)  and  (2.48)  followed  by  a simple  manipulation  of  (2.28),  (2.29),  and  (2.47) 
gives  the  following  first  order  differential  equations 
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(2.52) 


The  appropriate  boundary  conditions  for  the  above  equations  are 


(0)  = 1 

£(0  ) = /„ 


(2.31) 

(2.33) 


Tu(0)=  r“(0)-  r"(0)  = 1 

£■“(0)  = T^(O)-  r“(0)  = arbitrary  constant. 
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(2.53) 

(2.54) 


Here  the  two  variables  T“  and  E“  are  the  first  normal  stress  Vs  = r”  - and  the  hoop 
stress  E — — r” . The  above  system  (Eqs.(2.49)~(2.52))  can  be  integrated  by  using 

fourth  order  Runge-Kutta-Gill  method  along  with  the  boundary  conditions  (2.31),  (2.33), 
(2.53),  and  (2.54).  The  initial  values  of  the  hoop  stress  £"(0)  can  be  given  arbitrarily. 

The  results  of  the  computation  show  that  the  steady  state  solutions  for  the  axial 
velocity  and  the  stress  profiles  have  no  significant  difference  when  the  initial  value  of 
Ess( 0)  is  varied.  Figure  7 shows  the  hoop  stress  distributions  along  flow  direction  at  two 
different  skin  flow  rate  ratios  (f=0.2  and  0.5)  when  all  other  paramertes  are  specified  as 
Dr=  30,  De=  0.02,  M=  1,  /0=  1 0.  The  initial  hoop  stresses  are  varied  at  -10,  -2,  0,  2,  and 

10,  respectively.  In  spite  of  the  difference  in  the  initial  hoop  stress  condition  all  the  hoop 
stress  profiles  quickly  approach  to  zero  indicating  the  negligible  importance  of  hoop 
stress  in  fiber  spinning  flow.  Consequently,  the  system  we  are  dealing  with  may  not  be 

seriously  affected  by  the  initial  hoop  stress.  Thus,  the  initial  stress  condition 
Tw(°)  = Tr(°)  which  is  used  in  the  present  computation  should  not  seriously  diminish  the 

generality  of  the  current  study. 


2.6  Effect  of  Initial  Shear  Stress 

The  boundary  conditions  for  the  axial  velocity  described  in  section  2.3  are 
straightforward.  The  specification  of  the  initial  stress  state  of  the  viscoelastic  skin  layer, 
however,  is  somewhat  ambiguous  since  it  is  not  a measurable  quantity.  In  case  of  a 
purely  viscous  fluid,  the  initial  stress  state  does  not  depend  on  the  previous  strain  history 
and  the  specification  of  it  is  not  necessary.  The  stress  state  of  viscoelastic  fluids,  on  the 
other  hand,  depends  on  the  upstream  condition  due  to  their  memory  effect.  Therefore,  if 
the  stress  state  at  z=0  is  to  have  strong  influence  on  the  downstream  dynamics  of  the  fiber 


E(z)=Tee(z)-T„(z) 
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Figure  7.  Normal  stress  distribution  along  the  flow  direction  at 

different  initial  values  of  normal  stress:  (a)  f=  0.2,  Dr=  30, 
De= 0.02,  r0=  1 0,  M=  1;  (b)  f=  0.5,  Dr= 30,  0.02,  /0=10, 

M-l 


E(z)=T 9%(z)-T„(z) 
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Figure  7.  (continued). 
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spinning  flow,  it  will  be  necessary  to  consider  the  strain  history  of  the  material  inside  the 
die.  In  this  section,  influence  of  initial  stress  condition  (especially  to  value  of  /0)  on  the 
overall  dynamics  of  the  fiber  spinning  flow  is  investigated  by  integrating  the  system 
equations  for  various  values  of  t0. 

The  axial  normal  stress  distributions  for  two-phase  coextruded  fiber  (Newtonian 
and  UCM)  at  different  combinations  (f=  0.2,  0.5)  are  plotted  in  Figure  8a  and  8b, 

respectively.  It  is  apparent  from  the  figure  that  the  influence  of  the  initial  stress  is 
restricted  to  a rather  small  region  near  the  die  exit  (up  to  zsO.  1 ) and  it  has  only  negligible 

influence  in  most  of  the  draw-down  region.  This  behavior  is  conceivable  considering  the 
fact  that  the  fiber  spinning  flow  is  a "strong  flow"  in  which  the  deformation  of  the 
material  is  near  exponential.  Consequently,  within  a short  downstream  distance,  the  old 
strain  history  is  completely  overwhelmed  by  the  strong  downstream  conditions  (Larson, 
1988). 

The  velocity  distributions  given  in  Figure  9a  and  9b  also  reflect  the  negligible 
influence  of  the  initial  stress  condition.  Obviously,  the  influence  of  the  initial  stress  state 
will  persist  over  a longer  distance  if  t0  is  extremely  large. 


2.7  Effect  of  Shear  Viscosity  Ratio  M 

The  shear  viscosity  ratio  M is  one  of  the  parameters  which  affects  the  flow 
system.  For  the  two-phase  fiber  spinning  flow,  the  draw  ratio  as  a function  of  inverse 
axial  force  is  plotted  for  different  values  of  M (Figure  10).  It  may  be  easier  to  consider 
that  changing  the  value  of  M is  equivalent  to  changing  the  shear  viscosity  of  the 
Newtonian  core  layer  while  that  of  the  viscoelastic  skin  is  fixed.  That  is,  the  value  of  M 
is  decreased  as  the  shear  viscosity  of  the  Newtonian  core  decreases.  At  the  given 
conditions  (De= 0.02,  /0=  1 0),  the  Dr-G  curves  for  M=\  and/=0.2,  0.5  (curve  2 and  5 in 
Figure  10)  have  been  discussed  in  section  2.4. 


Figure  8.  Shear  stress  distribution  along  the  flow  direction  at  different 
initial  values  of  shear  stress:  (a)  f=  0.2,  Dr=  30,  De= 0.02, 
/0=10,  M=\\ (b)/=0.5,  Dr=  30,  De=  0.02,  /0=  1 0,  M=  1. 
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Figure  8.  (continued). 
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Figure  9.  Steady  state  axial  velocity  profiles  at  different  initial  values 
of  shear  stress:  (a)  f=  0.2,  Dr=  30,  De=  0.02,  /0=10,  M=  1; 
(b)  f=  0.5,  Dr=30,  De=  0.02,  t0=  10,  A/=l. 
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Figure  9.  (continue). 
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Figure  10.  Draw  ratio  as  a function  of  inverse  axial  force  at  different 

shear  viscosity  ratios;  De=0.02,  r0=10. 
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When  M is  small  (e.g.  A/=0.2,  curve  1)  the  relative  influence  of  the  Newtonian 
fluid  becomes  weaker  and  the  flow  dynamics  is  dictated  by  the  viscoelastic  skin  layer  at  a 
smaller  value  of  Dr  compared  to  larger  M cases  (i.e.,  curves  2 and  3).  Consequenly, 
curve  1 for  small  M does  not  show  the  multiple  steady  state  solutions  that  are  observed 
with  curve  2 or  3.  The  same  behavior  is  also  observed  with  a larger  values  of  / (i.e., 
curves  4,  5 and  6).  The  stability  behavior  of  these  curves  are  discussed  in  Chapter  3. 

2.8  Summary  and  Conclusion 

The  steady  state  flow  characteristics  of  a two-phase  fiber  spinning  flow  has  been 
investigated  in  this  Chapter.  It  has  been  assumed  that  the  core  and  the  skin  layers  are  a 
Newtonian  and  an  upper-convected  Maxwell  fluid,  respectively.  These  constitutive 
models  are  chosen  since  the  major  objective  of  the  study  is  to  investigate  the  competing 
influence  of  the  viscous  and  viscoelastic  forces  of  the  core  and  skin  layers.  In  order  to 
avoid  unnecessary  complexity  it  has  been  further  assumed  that  the  flow  is  isothermal  and 
that  all  external  and  surface  forces  are  negligible  other  than  the  viscous  and  viscoelastic 
forces.  Based  on  the  results  of  the  asymptotic  analyses  by  previous  investigators,  the  one- 
dimensional approximation  has  been  applied  at  the  outset  in  formulating  the  system 
equations.  The  one-dimensional  approximation  is  equivalent  to  assuming  that  the 
velocity  and  the  stress  fields  of  the  fiber  spinning  flow  is  independent  of  the  radial 
coordinate. 

Among  many  variables  which  may  affect  the  two-phase  flow,  the  effect  of  the 
flow  rate  ratio  f and  the  Deborah  number  De  has  been  extensively  studied.  It  is  rather 
obvious  to  anticipate  that  the  characteristics  of  the  two-phase  flow  may  sit  between  those 
of  the  single-phase  flows  of  each  constituent  fluid.  An  unexpected  behavior  has  been  also 
predicted  by  the  current  model  in  that  the  coextrusion  flow  may  enable  us  to  obtain  a 
higher  draw  ratio  than  possible  with  the  upper-convected  Maxwell  fluid  only.  As  it  is 
well  established  that  the  spinning  of  a Newtonian  fluid  is  limited  to  a relatively  low  value 
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of  draw  ratio  due  to  the  draw  resonance  instability,  this  new  finding  may  have  important 
implications  if  the  high  draw  ratio  flow  is  confirmed  to  be  stable.  In  the  following 

Chapter,  the  stability  of  the  two-phase  flow  is  investigated  by  the  method  of  linear 
stability  analysis. 

Although  the  system  under  current  consideration  is  very  simple,  the  procedure 
and  results  presented  here  may  serve  as  a basis  for  treating  more  complicated  elaborate 
two-phase  systems,  with  more  sophisticated  rheological  equation. 


CHAPTER  3 

LINEAR  STABILITY  ANALYSIS 


3.1  Introduction 

The  physical  mechanism  of  the  draw  resonance  instability  is  not  yet  clearly 
understood.  It,  however,  appears  to  be  a kinematic  instability  in  the  sense  that  it  is 
observed  and  predicted  with  purely  viscous  fluids  even  in  the  absence  of  all  external  and 
surface  forces. 

In  case  of  a Newtonian  fluid,  only  one  dimensionless  parameter  (draw  ratio  Dr) 
appears  in  the  system  equations  and  the  predicted  instability  is  independent  of  the  shear 
viscosity,  the  flow  rate  or  the  draw  span.  In  case  of  viscoelastic  fluids,  on  the  other  hand, 
the  draw  resonance  is  influenced  by  the  viscoelasticity  which  is  represented  by  the 
Deborah  number,  De.  Since  the  Deborah  number  is  the  ratio  of  the  material  time  constant 
to  the  processing  time,  the  change  in  either  the  flow  rate  or  the  draw  span  results  in  the 
variation  of  De  thus  affecting  the  onset  of  draw  resonance.  When  a purely  viscous  fluid 
and  a viscoelastic  fluid  are  spun  together  as  in  the  present  system,  a unique  situation  may 
arise  since  the  relative  influence  of  the  viscoelasticity  may  change  with  the  flow  rate  ratio 
of  the  two  materials  even  for  a fixed  value  of  De.  In  this  Chapter  the  stability  of  the  two- 

phase  spinning  flow  is  investigated  using  the  steady  state  solutions  given  in  Chapter  2 as 
the  base  flow. 

The  method  of  linear  stability  analysis  is  applied  to  the  present  system  for  which 
the  mathematical  treatment  is  straightforward  although  complicated.  Only  a limited 
information  may  be  obtained  by  the  linear  analysis  as  it  is  relevant  to  only  infinitesimal 
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disturbances.  The  sustained  periodic  nature  of  the  draw  resonance  instability,  however, 
may  make  the  linear  stability  analysis  meaningful.  While  the  accurate  prediction  of  the 
wavelength  and  the  amplitude  of  the  periodic  instability  may  be  obtained  only  through  a 
nonlinear  analysis,  the  mathematical  method  adopted  in  this  work  allows  us  to  estimate 
the  wavelength  of  the  draw  resonance.  The  estimated  trend  in  the  wavelength  of  the 
instability  is  apparently  in  qualitative  agreement  with  the  experimental  observations 
which  will  be  discussed  in  Chapter  4.  In  the  following  sections  the  existing  studies  on  the 
stability  of  single-phase  flows  are  briefly  reviewed  followed  by  the  detailed  formulation 
and  analysis  for  the  two-phase  system. 

3.2  Stability  of  Isothermal  Newtonian  Fiber  Spinning  Flow 

The  stability  analysis  for  the  single-phase  flow  of  a viscous  fluid  was  first 
performed  by  Kase  and  Matsuo  (1967),  Matovich  and  Pearson  (1969).  Their  analyses 
were  in  principle  a linear  stability  analysis  although  they  involved  some  mathematical 
ambiguities.  Later,  Gelder  (1971)  treated  the  same  problem  in  a more  comprehensive 
manner  by  using  the  finite  difference  method.  As  another  method  for  solving  the  stability 
problem,  the  orthogonal  collocation  method  was  applied  by  Gupta  and  Ballman  (1982). 
The  various  mathematical  approaches  predicted  that  the  critical  draw  ratio  for  the  onset 
of  draw  resonance  is  20.21  for  Newtonian  fluids.  That  is,  the  spinning  flow  of  a 
Newtonian  fluid  is  stable  up  to  Dr=20.2l.  Beyond  that  critical  value,  the  spinning  flow 
becomes  unstable. 

In  this  study,  the  finite  difference  method  of  Gelder  was  adopted  for  the  two- 
phase  flow  (see  Appendix  A).  This  method  involves  the  conversion  of  a differential 
eigenvalue  problem  into  an  algebraic  eigenvalue  problem  by  a finite  differencing 
scheme,  which  is  then  solved  by  a well  developed  standard  technique.  While  the 
conversion  of  the  problem  into  a matrix  form  is  rather  complicated  algebraically  and  very 
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tedious,  this  method  has  an  extra  advantage  of  being  able  to  investigate  higher  mode 
eigenvalues  in  addition  to  the  most  dangerous  one. 

As  a measure  of  accuracy  of  the  present  calculation,  the  single-phase  flow  of  a 
Newtonian  fluid  was  treated  and  its  results  compared  with  the  existing  one.  In  Table  1, 
the  present  computational  results  are  compared  with  those  of  Chang  et  al.  (1981)  for 
various  spatial  resolution.  While  the  slight  numerical  difference  may  result  from  the 
difference  in  the  differencing  scheme  (i.e.,  forward,  central,  backward  differencing),  this 
comparison  confirms  that  present  computational  procedure  is  correct  and  accurate. 

3.3  Stability  of  Isothermal  UCM  Fiber  Spinning  Flow 

For  viscoelastic  fluids,  the  stability  of  fiber  spinning  flow  was  first  studied  by 
Fisher  and  Denn  (1976).  The  generalized  upper-convected  Maxwell  model  (or  White- 
Metzner  model)  was  used  for  the  analysis  in  which  the  shear  viscosity  and  the  elastic 
modulus  are  allowed  to  be  dependent  on  the  shear  rate.  The  UCM  model  is  a special  case 
in  which  both  the  shear  viscosity  and  the  elastic  modulus  are  fixed  as  constants.  When 
f=  1,  the  present  two-phase  system  is  reduced  to  the  single-phase  flow  of  a UCM  fluid. 
Therefore,  the  present  results  can  be  also  compared  with  those  of  Fisher  and  Denn 
(1976)  for  this  extreme  case.  The  stability  analysis  shows  that  the  viscoelasticity 
generally  has  a stabilizing  influence  of  the  draw  resonance  especially  when  the  Deborah 
number  is  large. 


3.4  Stability  of  Two-Phase  Fiber  Spinning  Flow 
There  exists  no  study  on  the  stability  of  a two-phase  fiber  spinning  flow  to  date. 
As  mentioned  in  the  preceding  sections,  existing  analyses  are  only  for  the  single-phase 
flows.  The  stability  of  the  two-phase  fiber  spinning  system  is  investigated  in  this  section 
for  the  first  time. 


Table  1.  Values  of  A1  for  isothermal,  Newtonian  spinning  at  Dr=  20.21 


Discretisation  points 


imag 


N (Chang  et  al.,  1981)  (Present  calculation) 


11 

-1.0959 

12.594 

-1.1033 

12.3898 

21 

-0.1489 

13.561 

-0.1811 

13.5519 

31 

-0.0547 

13.808 

-0.0683 

13.8061 

41 

-0.0294 

13.895 

-0.0367 

13.8954 

51 

-0.0205 

13.930 

-0.0233 

13.9365 

N— ><*> 

0.0 

14.012 

0.0012 

14.0092 
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The  present  study  is  only  a linear  stability  analysis  and  the  results  indicate  that  the 
onset  of  draw  resonance  of  a Newtonian  fluid  is  delayed  to  a higher  draw  ratio  by  the 
presence  of  the  viscoelastic  skin  layer.  It  is  also  predicted  that  in  a coextrusion  flow 
under  certain  conditions,  the  highest  draw  ratio  for  a stable  spinning  can  be  much  larger 
than  that  of  the  single-phase  spinning  of  either  fluid  (Newtonian  or  UCM). 

The  stability  of  steady  state  solutions  to  infinitesimal  disturbances  is  determined 
by  the  following  standard  procedure  of  the  perturbation  analysis  in  which  solutions  of  the 
following  form  are  sought: 


A(z,t)  = Ass(z)[l  + a(z)eX/]  (3.1) 

w{z,t)  = wss(z)[\  + w'{z)eXj ] (3.2) 

tzz(z,t)  = tZ(z)[\+xa(z)eXl]  (3.3) 

T(z,t)  = Tss(z)[\  + T'(z)ex‘]  (3.4) 


Here  the  superscripts  ss  denotes  the  steady  state,  and  T is  defined  as  x - r . The 
perturbation  functions  a , w”,  and  T’  are  assumed  to  be  very  small,  and  the  real  part 

of  X,  which  is  the  eigenvalue  of  a normal  mode,  determines  the  stability  of  the  steady 
solution. 

By  substitutiting  (3. 1)— (3.4)  into  the  system  partial  differential  equations  (2.1), 
(2. 13), (2. 14)  and  (2.25)  and  by  linearizing  the  resulting  equations  about  the  perturbation 
functions,  the  following  homogeneous  ordinary  differential  equations  are  derived 


wa a ' +wss w' ' = -hi 


(3.5) 
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(3.7) 


(3.8) 


Here  the  prime  denotes  the  differentiation  with  respect  to  the  axial  coordinate  z.  For  the 
spinning  process  in  which  the  take-up  speed  is  held  constant,  the  boundary  conditions  for 
the  eigensystem  (i.e.,  Equations  (3.5)-(3.8))  follow  from  Equations  (2. 19)— (2.23): 

a (0)  = w*  (0)  = r;  (0)  = T (0)  = w*  (1)  = 0 (3.9) 


Equations  (3.5)— (3.9)  constitute  an  eigenvalue  problem  for  which  nontrivial  solutions 
exist  only  for  certain  values  of  A steady  state  solution  is  stable  to  infinitesimal 
disturbances  if  all  eigenvalues  \ have  negative  real  parts,  since  the  perturbations  defined 
by  (3. 1 ) — (3.4)  then  decay  to  zero. 
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To  discretize  the  eigenvalue  problem  (Equations  (3.5M3.9)),  the  finite  difference 
method  was  used  with  the  step  size  h.  The  grid  points  are  either  boundary  or  regular 
interior  points.  For  all  grid  points  in  the  interior,  the  central  difference  was  applied. 
Either  forward  difference  or  backward  difference  was  used  at  the  two  boundary  points. 
Using  this  finite  differencing  scheme,  the  system  equations  are  converted  into  a standard 
form  of  an  algebraic  eigenvalue  problem.  For  the  spatial  resolution  of  n,  the  dimension  of 
the  eigenvector  is  3 n.  The  detailed  derivation  of  the  algebraic  eigenvalue  problem  is 
described  in  Appendix  A.  The  eigenvalue  equations  are  solved  by  using  the  subroutine 
called  DEVLRG  in  IMSL  mathematical  package. 

3.5  Results  and  Discussion 

In  principle,  there  exist  infinite  number  of  eigenvalues  for  the  given  eigenvalue 
problem  indicating  that  the  number  of  discretisation  points  should  be  infinite.  Due  to  the 
computational  limitation,  the  calculation  can  be  done  only  with  a finite  number  of 
discretisation  points.  Then,  the  numerical  value  of  eigenvalues  will  depend  on  the  spatial 
resolution  N.  In  Table  2,  the  dependence  of  the  first  eigenvalue  (i.e.,  the  eigenvalue  with 
largest  real  part)  on  the  number  of  discretisation  points  N is  given  for /=0. 2 and  0.5  when 
Dr=30  and  De=  0.02.  The  Richardson  extrapolation  (Froberg,  1985)  is  applied  to  these 
data  to  obtain  the  limiting  value  of  X as  the  number  of  discretisation  points  tends  to 

infinity.  The  positive  real  part  for  f=0.2  indicates  that  the  steady  state  solution  at  Dr= 30 
is  unstable  whereas  it  is  stable  when  ^=0.5  since  the  most  dangerous  eigenvalue  has  the 
negative  real  part.  Thus,  in  Figure  6a,  the  two  velocity  profiles  for  /= 0.5  and  1.0  are 
stable  to  infinitesimal  disturbances  while  the  others  are  unstable. 

In  Figure  1 1 the  real  parts  of  the  first  eigenvalue  is  plotted  as  a function  of  draw 
ratio  Dr  for ^=0.5  at  three  different  values  of  Deborah  number.  This  figure  indicates  that 
when  De= 0.016  the  spinning  flow  is  stable  at  all  draw  ratios  smaller  than  the  maximum 
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Table  2.  Values  of  A1  for  differential  values  of  discretisation  points 

(De=0.02,  Dr=  30,  /o=10.0,  M=1.0) 


Discretisation  points 
N 

^real 

/=0.2 

y} 

^iniag 

/=0.5 

'Veal 

A' 

unag 

11 

-0.6521 

12.8830 

-0.9022 

11.5929 

21 

0.2415 

14.3838 

-0.9969 

14.0687 

31 

0.2959 

14.7408 

-1.2054 

14.6858 

41 

0.2938 

14.8688 

-1.3359 

14.8928 

51 

0.2977 

14.9406 

-1.3841 

15.0186 

61 

0.2992 

14.9879 

-1.4117 

15.1019 

71 

0.2994 

15.0209 

-1.4208 

15.1591 

N— >oo 

0.3211 

15.1006 

-1.4272 

15.3224 
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Figure  11.  Variation  of  the  real  part  of  the  first  eigenvalue  with  draw 

ratio  for  various  Deborah  numbers  at  f=  0.5,  /0=10,  M=\ 
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attainable  value  which  corresponds  to  the  infinite  applied  force.  When  0.012  or 
0.014,  on  the  other  hand,  there  exist  two  stable  regions:  one  at  a low  value  of  Dr  and  the 
other  at  a higher  value  of  Dr.  These  behaviors  are  very  similar  to  those  of  a single-phase 
flow  of  an  UCM  fluid  for  which  the  upper  stable  region  at  high  draw  ratio  was  first 
found  by  Fisher  and  Denn  (1976).  The  general  behavior  given  in  Figure  1 1 is  typical  of 
all /in  the  range  of  between  0.1  to  0.8. 

The  critical  draw  ratio  at  which  the  real  part  of  the  first  eigenvalue  becomes  zero 
is  plotted  as  a function  of  Deborah  number  in  Figure  12.  In  this  plot  of  Dr*  vs.  De,  each 
pair  of  neutrally  stable  points  described  in  Figure  1 1 appears  as  a continuous  line 
dividing  the  whole  domain  into  the  stable  and  unstable  regions.  In  typical  fiber  spinning 
operations,  the  spinning  conditions  (i.e.,  output,  draw  span,  material,  etc.)  are 
predetermined  and  the  take-up  speed  is  gradually  increased  to  achieve  a desired  draw 
ratio.  In  other  word,  Dr  is  gradually  increased  to  the  target  value  for  a fixed  De. 
Therefore,  the  upper  stable  region  is  generally  unattainable  unless  special  approach  is 
taken  to  go  around  the  nose  (or  tip)  of  the  neutral-stability  curve.  While  a few 
investigators  alleges  experimental  observation  of  the  upper  stable  region  (Fisher  and 
Denn,  1976),  others  contend  that  it  is  an  artifact  of  the  unrealistic  constitutive  model.  In 
any  case,  the  more  important  character  of  the  neutral-stability  curve  is  the  sharp  upward 
turn  from  a near  horizontal  shape  with  increasing  De  whether  it  turns  further  backward 
forming  the  upper  stable  region  or  not.  This  sharp  upward  turn  occurs  at  about  £>e=0.01 
when  j=  1.0.  When  Dr  reaches  about  100,  the  fiber  may  snap  instead.  Thus,  when  De> 0.1, 
the  draw  resonance  is  unobservable  (or  the  critical  draw  ratio  for  its  onset  increases 
drastically).  The  dotted  diagonal  line  (c)  indicates  the  maximum  obtainable  draw  ratio  for 
fr  10  at  which  the  axial  tension  becomes  infinity.  This  behavior  is  physically  plausible 
since  the  viscoelasticity  effect  is  influential  only  when  the  Deborah  number  is  large. 

As  the  skin-layer  flow  rate  fraction  / is  decreased,  the  nose  of  the  unstable 
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Figure  12.  Critical  draw  ratio  as  a function  of  Deborah  number  for 

various  / at  f0=10,  M=\  (maximum  attainable  Dr  for  (a) 
/=0.2(Max);  (b)/=0.2(Asymptote);  ( c)J=0.5  and  1.0) 
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envelope  extends  to  a higher  value  of  De.  When  De  is  small,  the  lower  critical  draw  ratio 
does  not  vary  significantly  with  / and  it  is  close  to  the  Newtonian  value  of  20.21.  At  a 
higher  De  (i.e.,  De>0.01),  however,  the  lower  critical  draw  ratio  increases  with  / 
indicating  the  stabilizing  effect  of  the  viscoelastic  skin  layer.  From  this  figure  the 
stability  of  the  velocity  profiles  given  in  Figure  6 can  be  determined.  The  steady  state 
solution  described  in  Figure  6a  (De=0.02  and  Dr=3Q)  falls  on  the  stable  region  if^=1.0  or 
0.5  whereas  it  falls  on  the  unstable  region  ifj=0.2.  When  Dr  is  increased  up  to  60  at  the 
same  De  (i.e.,  De= 0.02),  the  spinning  condition  belongs  to  the  upper  stable  region  of 
/=0.2  curve.  Thus,  in  Figure  6 the  two  curves  for ^1.0  and  0.5  are  stable  whereas  those 
for  f=  0.2  and  0.0  are  not.  It  should  be  pointed  out  again  that  the  upper  stable  region  may 
be  an  artifact  which  may  not  exist  in  real  systems. 

The  most  interesting  feature  of  the  coextrusion  fiber  spinning  may  be  the 
possibility  of  a substantial  increase  in  the  maximum  attainable  draw  ratio  when  / is 
relatively  small.  As  discussed  previously,  Figure  13  indicates  that  a draw  ratio  which  is 
much  greater  than  the  asymptotic  value  for  an  infinite  draw  force  (i.e.,  G->0)  is 

attainable  when^=0.2.  The  stability  analysis  shows  that  the  low-G  branch  of  the  multiple 
steady  solutions  (i.e.,  the  portion  of  the/=0.2  curve  where  G<0.045  including  point  a)  is 
always  unstable.  The  high-G  branch  of  the  solution,  on  the  other  hand,  is  either  stable  or 
goes  through  two  transitions  with  increasing  Dr  depending  upon  the  magnitude  of  De. 
These  behaviors  are  clearly  shown  in  Figure  12  and  13.  When  De= 0.02  and  f= 0.2,  a 
stable-to-unstable  and  an  unstable-to-stable  transitions  occur  at  Dr=24.07  and  51.24, 
respectively  as  Dr  is  increased.  When  De  is  greater  than  0.0236,  on  the  other  hand,  the 
spinning  condition  is  outside  the  nose  of  the  unstable  envelope  and  the  "spinning  stability 
is  sustained  up  to  62.  The  solid  line  (a)  in  Figure  12  is  the  locus  of  the  maximum  points 
of  Dr  vs.  G curve  (e.g.  point  c in  Figure  13)  for  f=  0.2  which  includes 
(De,Dr)=( 0.0236,62)  point.  When/=1.0  or  0.5,  the  maximum  attainable  draw  ratio  at 
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Figure  13.  Draw  ratio  as  a function  of  inverse  axial  force  (/= 0.2, 

De=0.02,  /0=10,  M=  1) 
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De=0.0236  is  about  43.  Therefore,  much  higher  value  of  Dr  is  obtainable  when  f=  0.2 
than  when /0.5  or  1 .0. 

In  Figure  14,  Dr  vs.  G curves  are  given  for  various  values  of  De  when  /=0.2.  Also 
indicated  in  the  figure  are  the  neutral-stability  curves  dividing  the  stable  and  the  unstable 
region.  It  may  be  noted  that  the  neutral-stability  curve  on  the  left  hand  side  coincides 
with  the  locus  of  the  maximum  points  of  Dr-G  curves.  It  is  interesting  to  point  out  that 

the  imaginary  part  of  the  first  eigenvalue  is  zero  along  this  curve  (see  point  c in  Figure 
13)  indicating  catastrophic  failure  (e.g.,  sudden  snap  or  breakage)  beyond  this  curve.  The 
location  of  this  curve  depends  on  the  value  of / As /increases,  the  maximum  point  of  the 
Dr-G  curve  disappears  and  the  neutral  stability  curve  coincides  with  the  vertical  axia 

(i.e.,  Dr  axis).  On  the  other  hand,  the  neutral  stability  curve  on  the  right  hand  side  is  a 
universal  curve  whose  location  is  fixed  regardless  of  other  parameters.  This  right  hand 
side  neutral  stability  curve  cuts  through  the  portion  of  Dr  vs.  G curves  for  De< 0.0236. 
Since  in  a typical  operation,  the  high  Dr  is  approached  along  the  Dr-G  curve  from  the 

large-G  side,  the  stable  region  between  the  two  unstable  regions  is  not  attainable  when 
£>e<0.0236.  Therefore,  only  the  shaded  region  is  the  stable  region  that  is  attainable. 

One  of  the  most  important  findings  may  now  be  reiterated  as  follows.  For  a given 
choice  of  material,  the  Deborah  number  is  the  same  regardless  of  the  skin-layer  fraction/ 
if  the  total  flow  rate  is  the  same  and  the  spinning  occurs  over  the  same  draw  span.  Thus, 

it  is  possible  to  obtain  a much  higher  draw  ratio  with  coextrusion  than  obtainable  with 
either  the  UCM  fluid  or  the  Newtonian  fluid,  if  De  is  large  enough  so  that  the  high-G 

branch  solution  is  stable  without  cutting  through  the  unstable  region.  For  example,  if  a 
skin-layer  material  is  chosen  so  that  De  is  about  0.0236,  a stable  spinning  is  possible  up 
to  Dr=62  when/=0.2.  This  draw  ratio  is  much  larger  than  the  value  which  is  attainable 
with  either  the  core  material  (Drmax=20.21)  or  the  skin  material  (Drmax=43.37)  for  the 


same  value  of  De. 
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Figure  14.  Draw  ratio  as  a function  of  inverse  axial  force  for  various 
Deborah  numbers  (^=0.2,  /0=10,  M=\) 
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These  stability  characteristics  of  the  two-phase  flow  may  be  presented  in  another 
way  using  the  7>G\r\Dr  vs.  De  curves  discussed  in  Chapter  2.  These  have  been  already 
indicated  in  Figures  4a  and  4b.  When  f=\.0,  0.5  or  0.2,  the  system  for  Dr=  20  is  always 
stable  for  any  De  which  is  smaller  than  the  maximum  values.  It  is  not  surprising  when 
considering  the  critical  draw  ratio  of  20.21  for  a Newtonian  fluid.  An  interesting  point 
though  is  the  value  of  the  maximum  De  which  is  0.061  for f=0.2  that  is  much  larger  than 
those  for  /=0.5  or  1.0.  When  Dr  is  larger  (i.e.,  Figure  4b  for  Dr=  30),  the  system  is 
unstable  for  small  De  and  moves  into  a stable  region  as  De  is  increased.  Following  the 
3 GlnDr  vs.  De  curve  from  left  to  right  is  equivalent  to  following  a horizontal  line  in 
Figure  12  (i.e.,  constant  Dr  line).  Thus,  the  lower  neutral-stability  curve  in  Figure  15 
which  connects  the  maximum  De  points  is  equivalent  to  the  diagonal  line  (a)  in  Figure 
12.  The  unstable  region  in  the  upper  left  comer  is  then  the  draw  resonance  region  and  the 
lower  unstable  region  is  the  catastrophic  failure  (i.e.,  breakage)  region. 

An  important  advantage  of  using  the  finite  difference  technique  for  the 
eigenvalue  problem  is  the  information  regarding  the  higher  modes  of  instability  (Geyling 
and  Homsy,  1980).  In  Table  3 the  real  parts  of  the  first  four  eigenvalues  are  given  at  the 
two  critical  draw  ratios.  These  data  represent  the  general  trend  of  eigenvalues  at  all 
critical  conditions  that  have  been  examined.  It  may  be  noted  that  the  higher  modes  for 
these  cases  have  relatively  large  negative  values  (i.e.,  Areal<-2.6545)  indicating  that  they 

are  strongly  damped.  Although  the  linear  stability  analysis  provides  only  a limited 
information  about  the  dynamical  behavior  of  the  process  for  small  time,  this  trend  in  the 
real  part  of  the  eigenvalues  suggests  that  the  instability  may  quickly  evolve  to  a nonlinear 
regime  with  the  wavelength  (or  the  period)  defined  by  the  imaginary  part  of  the  first 
eigenvalue  (i.e.,  14.05  for  Dr=21.79  and  17.78  for  Dr=75.96).  Thus,  a measured  period 
of  draw  resonance  instability  may  agree  with  this  predicted  value  if  the  constitutive 
model  properly  represents  the  spinning  material. 


3G(lnDr) 
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Figure  15.  Inverse  axial  force  as  a function  of  Deborah  number  for 

various  draw  ratios  (^=0.2,  /o=10,  M=\) 


Table  3.  Real  and  imaginary  parts  of  first  four  eigenvalues  at  the 
critical  Dr * (f=  0.2,  De=0.015,  /0=10.0,  A/=1.0) 


Dr *=21. 79 


Dr *=75.96 


^"real 

Ximag 

\«al 

Ximag 

1st  eigenvalue 

0.0000 

14.0517 

0.000 

17.7820 

2nd  eigenvalue 

-2.6545 

33.8660 

-7.1597 

47.2597 

3rd  eigenvalue 

-4.5370 

62.9363 

-10.6844 

67.6269 

4th  eigenvalue 

-4.7455 

49.9121 

-14.2681 

83.1421 
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In  Figure  16  the  neutral-stability  curve  for/=0.2  is  given  along  which  all  the  real 

parts  of  the  first  eigenvalues  are  zero.  Also  given  on  the  curve  are  the  imaginary  parts  of 
the  first  eigenvalues,  A'imag.  A noticeable  trend  is  that  the  imaginary  part  of  the  first 

eigenvalue  increases  with  the  critical  draw  ratio,  Dr*. 

The  dimensionless  period  t (scaled  by  L/w0)  which  corresponds  to  A1  can  be 

If 

calculated  by  the  following  equation. 

2k 

t = ir-  (3.10) 

imag 

In  Figure  17,  the  period  thus  calculated  is  given  as  a function  of  critical  draw  ratio  for 

various  values  of  f.  This  figure  indicates  that  the  period  (thus  the  wavelength  of  the 

instability)  decreases  with  increasing  critical  draw  ratio.  As  pointed  out  earlier  this 

predicted  period  (or  wavelength)  may  agree  with  measured  value  if  the  constitutive 

model  accurately  describes  the  rheology  of  the  material  being  spun.  The  UCM  fluid 

model  is  too  simple  to  be  realistic.  Thus,  quantitative  agreement  between  the  measured 
and  predicted  period  is  not  expected.  The  Newtonian  limit  (i.e.,  A1.  =14),  however, 

should  show  the  agreement  if  the  current  argument  is  correct.  Unfortunately,  no 
polymers  are  Newtonian.  However,  the  relative  insensitivity  of  the  critical  Dr*  over  a 
broad  range  of  De  (see  Figure  16  for  £>e<0.01)  suggests  that  if  a polymer  with  a small 
relaxation  time  is  spun  the  period  of  draw  resonance  instability  may  closely  simulate  the 
Newtonain  value.  As  will  be  discussed  in  Chapter  4,  the  quantitative  agreement  is  in  fact 
obtained  in  the  Newtonian  limit.  It  will  be  also  shown  that  the  decreasing  trend  of  period 
with  increasing  Dr*  is  also  in  accordance  with  the  experimental  observations. 

A characteristic  time  which  may  be  closely  related  with  the  kinematics  of  the 
fiber  spinning  flow  is  the  material  residence  time,  7^r.  Trst  is  the  time  needed  for  a 
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Figure  16.  Neutral-stability  curve  for/=0.2,  /0=10,  M=  1.  Values  of 

A1-  . indicated  on  curve. 


Period 
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Figure  17.  Period  of  the  first  normal  mode  at  critical  conditions. 


Values  of  Alimjg  indicated  on  curve. 
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small  fluid  element  to  travel  from  the  die  exit  to  the  take-up  (i.e.,  Lagrangian  time  at  the 
take-up).  Once  the  axial  velocity  is  determined,  the  material  residence  time  can  be 
calculated  as  follows: 


As  discussed  in  Chapter  2,  the  axial  velocity  profile  is  significantly  affected  by  the 
viscoelasticity  through  the  changes  in  the  kinematics  of  the  flow.  If  the  draw  resonance 
instability  is  to  be  closely  related  to  the  kinematics  of  the  flow,  so  is  the  characteristic 
time  of  the  instability  (i.e.,  period)  to  the  residence  time  which  may  represent  the  flow 
kinematics.  In  Figure  18,  the  dimensionless  period  t at  the  critical  draw  ratio  is  plotted 
against  the  residence  time  determined  from  the  steady  velocity  profile  at  that  draw  ratio. 
The  Simpson's  3/8  integration  rule  has  been  used  to  calculate  the  residence  time. 
Although  it  is  not  clear  at  this  point  what  the  relationship  between  t and  Trst  should  be 
if  there  exists  any,  Figure  18  shows  a rather  strong  correlation  between  the  two 
characteristic  times  in  that  the  period  of  the  instability  is  small  if  the  material  residence 
time  is  short.  All  data  are  near  the  straight  dashed  line  which  is  a linear  regression  from 

those  data  of/=1.0  (i.e.,  UCM  single-phase  fiber).  The  regression  expression  is  written 
as  follows: 


/ = 0.062+1. 226  xTkst 


(3.12) 


The  influence  of  the  shear  viscosity  ratio  M on  the  steady  spinning  flow  has  been 
discussed  in  the  previous  Chapter.  Here  the  influence  of  M on  the  flow  stability  has  been 
investigated  and  the  results  are  given  in  Figure  19.  For  a fixed  value  of / (i.e.,/=0.2),  the 
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Figure  18.  Dimensionless  period  as  a function  of  residence  time  for 

various / at  /0=10  and  M=\ 


68 


De 


Figure  19.  Critical  draw  ratio  as  a function  of  Deborah  number  for 

various  M at  /=0.2,  /0=10,  compared  with  that  for 
/=1.0(UCM),  M=\  and  /0=10;  curves:  neutral-stability 
curves,  straight  lines:  maximum  draw  ratio  locus 
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shear  viscosity  of  the  Newtonian  core  was  varies  to  get  M=0.2,  1.0  and  5.0,  respectively. 
The  neutral-stability  curves  shown  in  Figure  19  with  various  M show  that  when  the 
values  of  M decrease,  the  nose  of  the  neutral-stability  curve  on  the  Dr*-De  plane  is 

moved  toward  left.  This  trend  is  plausible  since  the  influence  of  the  Newtonian  core 
weakens  as  M decreases  resulting  in  the  shift  of  the  neutral-stability  curve  toward  that  for 
the  UCM  fluid.  Decreasing  M is  equivalent  to  decreasing  the  shear  viscosity  of  the 
Newtonian  core.  While  the  overall  shape  (or  trend)  of  the  neutral-stability  curve  remains 
the  same,  this  figure  bears  an  interesting  implication.  For  a fixed  total  flow  rate  and / if 
the  same  skin-layer  material  is  used  while  changing  the  viscosity  of  the  core  material,  the 
Deborah  number  remains  the  same.  If  De  is  0.02,  the  coextrusion  flow  for  M= 0.2  is 
stable  until  it  breaks  at  Dt=  50.  Draw  resonance,  however,  is  encountered  at  a much 

smaller  Dr  for  others  (i.e.,  M=  1.0  or  5.0).  If  the  skin  material  has  a large  relaxation  time 
so  that  De=0.4,  on  the  other  hand,  all  three  curves  indicate  stability.  The  maximum 

attainable  Dr  for  A/=5.0,  however,  is  about  60  whereas  those  for  others  (i.e.,  M- 0.2  or 
1.0)  are  smaller  than  about  35.  Thus,  from  a processing  standpoint  this  plot  suggests  that 
a higher  viscosity  of  the  Newtonian  core  is  favorable  as  long  as  the  skin  layer  has  a large 
enough  relaxation  time. 


3.6  Summary  and  Conclusion 

The  stability  of  a two-phase  fiber  spinning  flow  to  an  infinitesimal  disturbance 
has  been  investigated  in  this  Chapter.  By  the  method  of  linear  stability  analysis,  the 
critical  draw  ratio  for  the  onset  of  draw  resonance  has  been  determined  as  a function  of 
Deborah  number  De  and  the  skin-layer  flow  rate  fraction  f The  shape  of  the  neutral- 
stability  curve  in  a Dr  vs.  De  plane  for  the  two-phase  flow  is  similar  to  that  for  the 
single-phase  flow  of  a UCM  fluid. 
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When  De  is  small  (i.e.,  when  the  viscoelasticity  of  the  skin  layer  is  small),  the 
critical  draw  ratio  remains  in  the  vicinity  of  20.21  which  is  the  value  for  a Newtonian 
fluid.  As  De  is  increased,  the  neutral-stability  curve  make  a sharp  upward  turn  thus 
making  the  spinning  flow  free  from  the  draw  resonance  instability.  In  this  case,  the 
spinning  flow  is  limited  by  the  large  axial  stress  which  can  exceed  the  cohesive  strength 
of  the  material  resulting  in  a spinline  breakage.  This  general  trend  is  preserved  for  all 
values  of/although  the  value  of  De  where  the  sharp  upturn  occurs  shifts  to  a larger  value 
with  decreasing/  On  the  other  hand,  the  upper  limit  of  draw  ratio  at  which  the  spinline 
may  break  increases  with  decreasing  f Consequently,  an  interesting  behavior  can  be 
observed  with  the  coextrusion  flow  in  that  the  highest  draw  ratio  for  a stable  two-phase 
spinning  can  be  larger  than  that  for  each  constituent  material  under  a certain  condition. 
This  prediction  may  have  important  practical  implication  if  it  can  be  proved 
experimentally.  Although  the  present  result  is  based  on  a simple  constitutive  model,  the 
overall  characteristics  is  physically  plausible  and  other  viscoelastic  constitutive  model 
may  result  in  a similar  prediction. 

The  linear  stability  analysis  provides  only  a limited  information  for  the  short-time 
behavior  of  the  instability.  Nevertheless,  the  qualitative  trend  predicted  by  the  theory 
regarding  the  critical  draw  ratio  as  well  as  the  wavelength  of  the  instability  is  apparently 

in  reasonable  agreement  with  the  experimental  observations  which  will  be  discussed  in 
the  following  Chapter. 


CHAPTER  4 

DRAW  RESONANCE  EXPERIMENTS 

4.1  Introduction 

Coextrusion  fiber  spinning  experiments  have  been  conducted  using  a linear  low- 
density  polyethylene  (LLDPE)  and  a low-density  polyethylene  (LDPE)  for  the  core  and 
the  skin  layers,  respectively.  These  two  materials  are  chemically  identical,  yet  their 
rheological  properties  are  very  different  from  one  another  as  they  are  manufactured  by 
completely  different  processes.  The  LDPE  has  a high  degree  of  long  chain  branching, 
and  has  a large  relaxation  time  compared  to  the  LLDPE.  Consequently,  the  former  tends 
to  show  a distinct  strain  hardening  behavior  when  subjected  to  an  extensional  flow 
whereas  the  latter  does  not.  Furthermore,  the  LLDPE  is  known  to  be  very  prone  to  draw 
resonance  whereas  the  LDPE  is  not.  The  spinability  of  the  LDPE  is  in  fact  observed  to  be 
limited  by  the  breakage  of  the  spinline  rather  than  by  the  draw  resonance. 

The  objective  of  the  experimental  study  is  to  investigate  the  draw  resonance 
characteristics  of  a bicomponent  fiber  spinning  process  and  to  compare  them  with  the 
theoretical  investigation  described  in  the  previous  chapters. 

The  theory  mentioned  in  previous  chapters  is  for  a simple  constitutive  model  with 
many  simplifying  assumptions.  Consequently,  any  quantitative  agreement  with 
experimental  observations  may  not  be  expected.  Nevertheless,  the  general  trend  predicted 

by  the  theory  seems  to  be  in  reasonable  agreement  with  the  experimental  results  at  least 
qualitatively. 
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In  most  of  the  experiments  on  the  single-phase  fiber  spinning,  the  ratio  of 
maximum-to-minimum  diameter,  critical  draw  ratio  and  period  (and  sometimes  draw 
force)  were  measured  while  the  extrusion  output  (or  flow  rate),  take-up  velocity, 
spinning  length  (or  draw  span)  and  melt  temperature  were  varied.  Although  there  are 
some  contradicting  reports  depending  upon  the  choice  of  the  materials  and  extrusion 
conditions,  many  experiments  have  shown  the  following  general  characteristics  of  draw 
resonance: 

(1)  A longer  draw  span  results  in  a decrease  in  the  maximum-to-minimum  diameter  ratio 

and  an  increase  in  period,  whereas  it  has  negligible  influence  on  the  critical  draw 
ratio. 

(2)  Increasing  output  typically  has  the  opposite  effect  of  increasing  draw  span. 

(3)  A higher  melt  temperature  results  in  a larger  critical  draw  ratio  and  a smaller 
maximum-to-minimum  diameter  ratio. 

The  proper  time  scale  for  the  fiber  spinning  flow  is  the  processing  characteristic 
time  L/w0  which  is  the  ratio  of  the  draw  span  to  the  initial  average  velocity  as  discussed 
in  the  previous  chapters.  It  seems  that  the  observations  (1)  and  (2)  simply  reflect  the 
nature  of  the  time  scale  for  this  process. 

A similar  experimental  procedure  was  followed  for  the  present  study  in  which  the 
critical  draw  ratio,  the  ratio  of  the  maximum-to-minimum  diameter,  and  the  period  were 

measured  while  the  spinning  length  (or  draw  span)  and  the  flow  rate  fraction  of  the  skin 
layer  were  varied. 

4.2  Experiment  Materials 

Four  different  polyethylenes  were  used  for  the  experimental  study.  Two  of  them 
were  linear  low-density  polyethylene  (LLDPE)  and  the  rest  were  low-density 
polyethylene  (LDPE).  The  two  LLDPEs  were  for  the  core  while  the  LDPEs  were  for  the 
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skin  layer.  These  materials  are  commercial  products  of  Union  Carbide  Corporation  and 
DuPont.  As  pointed  out  previously,  LLDPEs  are  very  prone  to  draw  resonance  whereas 
LDPEs  are  not  due  to  their  large  relaxation  time.  The  significant  difference  in  their 
rheological  properties  results  from  their  different  molecular  structures.  None  of  these 
materials  may  be  adequately  described  by  either  the  Newtonian  or  the  upper-convected 
Maxwell  fluid  model.  Nevertheless,  these  materials  were  chosen  with  a Newtonian  and  a 
UCM  fluid  in  mind.  Despite  the  differences  in  their  rheological  properties,  all  of  them 
are  chemically  compatible  with  each  other  and  consequently,  good  adhesion  at  the  core- 
skin interface  is  ensured. 

In  Table  4,  the  nominal  information  for  the  materials  is  given  which  has  been 
provided  by  the  Manufacturers.  Also  included  are  the  zero-shear  viscosities  that  have 
been  measured  in  our  Department.  The  important  difference  between  the  two  LLDPEs  is 
their  shear  viscosities  represented  by  the  melt  index.  The  melt  index,  which  is  one  of  the 
standard  parameters  that  characterize  the  material,  represents  the  flow  rate  of  the  material 
through  a capillary  at  a fixed  temperature  and  a fixed  pressure  drop.  Thus,  a high  value 
of  melt  index  indicates  a low  viscosity.  In  case  of  LLDPEs,  the  shear  viscosity  is,  very 
roughly,  inversely  proportional  to  the  melt  index,  and  the  zero-shear  viscosities  seem  to 
confirm  the  rough  relationship. 

The  two  LDPEs  listed  in  Table  4 are  known  to  have  different  level  of  long  chain 
branching.  DYNH-1  (Union  Carbide)  is  produced  by  a tubular  reator  whereas  Alathon- 
1540  (DuPont)  is  by  an  autoclave  reactor.  Autoclave  resins  are  known  to  have  much 
higher  level  of  long  chain  branching  than  others,  and  consequently  have  larger  relaxation 
time  than  others  at  a similar  molecular  weight.  Thus,  it  is  expected  that  Alathon- 1 540 
will  have  a higher  Deborah  number  than  DYNH-1  under  the  same  extrusion  condition 
despite  its  larger  melt  index. 
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Table  4.  The  properties  of  experiment  materials 


Material 

type 

Nomenclature 

(Manufacturer) 

Melt 

Index 

Density,  p 
(kg/m3) 

Zero  shear 
viscosity, T)0 

(Pa-sec) 

(at  °C) 

LLDPE 

HS-7028 
(Union  Carbide) 

1 

918 

9919 

(200) 

LLDPE 

DMDA-8320 
(Union  Carbide) 

20 

922 

2228 

(200) 

LDPE 

DYNH-1 
(Union  Carbide) 

1 

918 

22240 

(145) 

LDPE 

Alathon-1540 

(DuPont) 

6.5 

918 

7899 

(145) 
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The  most  desirable  choice  of  materials  for  the  current  experiment  may  be  two 
LLDPEs  with  different  shear  viscosity  with  all  other  proporties  the  same  and  two  LDPEs 
with  different  elasticity  while  all  other  properties  remain  the  same.  In  case  of  LLDPEs, 
such  condition  is  met  reasonable  well  with  the  current  materials.  In  case  of  LDPEs, 
however,  the  shear  viscosities  are  very  different.  It  is  an  unfortunate  situation,  but  limited 
availability  of  materials  provided  no  other  choice. 

The  rheological  property  information  given  above  is  rather  qualitative.  The  most 
important  rheological  property  relevant  to  the  fiber  spinning  process  is  the  elongational 
viscosity.  The  measurement  of  steady  elongational  viscosity,  however,  is  very  difficult  if 
not  impossible.  Furthermore,  the  interpretation  of  the  data  provided  by  any  existing 
elongational  viscometer"  is  controversial  since  what  is  measured  may  not  be  a steady 
behavior  but  the  transient  behavior  of  the  strong  flow. 

In  Figures  20-23,  additional  information  is  given  for  the  rheological  properties  of 
the  materials.  These  were  measured  using  a rheogoniometer  (Model  # RMS-800). 
Parallel  plates  of  25  mm  in  diameter  were  ran  at  a steady  mode  to  measure  the  shear 
viscosity  and  the  normal  stress  differences.  The  gap  distance  were  varied  at  0.4,  0.8,  1.0 
1.1  and  1.5  mm.  The  temperature  for  the  measurements  was  fixed  at  145°C  and  200°C 
for  the  LDPEs  and  the  LLDPEs,  respectively,  to  match  with  the  operation  temperatures 
for  the  fiber  spinning  experiment. 

The  results  in  Figure  20  show  that  the  viscosities  of  all  materials  have  no  apparent 
change  with  shear  rate  at  low  shear  rate  range  whereas  DYNH-1  shows  a little  shear 
thinning  behavior  at  a higher  shear  rate.  The  high  viscosities  of  LDPEs  which  are 
apparently  higher  or  equivalent  to  that  of  HS-7028  are  due  to  the  low  measurement 
temperature.  The  DMDA-8320  has  a much  smaller  viscosity  than  others  as  expected. 
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Figure  20.  Shear  viscosity  function  r|  versus  shear  rate  y of  linear 

low-density  polyethylene  at  200°C:  (o)  HS-7028;  (□) 

DMDA-8320,  and  low-density  polyethylene  at  145°C:  (a) 
DYNH-1;  (*)  Alathon-1540. 
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Shear  rate  (sec-1) 


Figure  21.  Shear  stress  x12  versus  shear  rate  y of  linear  low-density 

polyethylene  at  200°C:  (o)  HS-7028;  (□)  DMDA-8320, 
and  low-density  polyethylene  at  145°C:  (a)  DYNH-1;  (*) 
Alathon-1540. 


Normal  stress  difference 
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Figure  22.  Normal  stress  difference  Nl-N2  versus  shear  rate  y of  linear 

low-density  polyethylene  at  200°C:  (o)  HS-7028;  (□) 
DMDA-8320,  and  low-density  polyethylene  at  145°C:  (a) 
DYNH-1;  (*)  Alathon-1540. 


Normal  stress  function. *!-*,  (Pa-sec*) 
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Figure  23.  Normal  stress  function  xVl-y¥2=(Nl-N2)/y2  versus  shear 

rate  y of  linear  low-density  polyethylene  at  200°C:  (o) 

HS-7028;  (□)  DMDA-8320,  and  low-density  polyethylene 
at  145°C:  (a)  DYNH-1;  (*)  Alathon-1540. 
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The  difference  between  the  first  and  the  second  normal  stress  difference  is  plotted 
as  a function  of  shear  rate  in  Figure  22  and  Figure  23  represents  the  difference  between 
the  normal  stress  functions,  'P,-'P2.  The  normal  stress  function  MC  is  defined  as 
= iV.  / y2  0=1,2).  The  normal  pressure  on  the  plate  at  the  center  or  at  the  rim  cannot  be 
measured  with  the  current  instrument.  Thus,  the  first  and  the  second  normal  stresses 
could  not  be  decoupled.  The  second  normal  stress  difference  N2,  however,  is  typically 
much  smaller  than  the  first  normal  stress  difference,  and  (NrN2)  may  well  represent  Nv 

It  was  pointed  out  that  Alathon-1540  might  have  higher  viscoelasticity  than  DYNH-1  due 
to  its  higher  degree  of  long  chain  branching.  Contrary  to  the  speculation,  (N,-A2)  of 

DYNH-1  is  larger  than  that  of  Alathon-1540.  It  may  be  due  to  the  high  shear  viscosity  of 
DYNH-1.  The  first  normal  stress  function  VF,(=  2V,  / y2)  of  a UCM  fluid  is  proportional 
to  its  shear  viscosity;  that  is 

'F,  = 2rjA  (4.1) 

Although  none  of  the  LDPEs  may  be  well  represented  by  the  UCM  model,  this  equation 
indicates  the  strong  influence  of  shear  viscosity  on  T*,. 

All  these  properties  are  for  a shear  flow  situation  and  the  direct  correlation 
between  the  shear  flow  properties  and  those  for  an  elongational  flow  is  unknown. 
However,  if  it  is  assumed  that  the  first  normal  stress  difference  may  represent  the 
magnitude  of  the  material  time  constant  (or  relaxation  time)  X,  DYNH-1  may  have  a 
larger  relaxation  time  than  Alathon- 1 540. 

Even  though  the  shear  viscosity  of  HS-7028  is  much  larger  than  Alathon-1540,  its 
(Nl-N2)  is  much  smaller  indicating  the  low  viscoelasticity  of  HS-7028.  As  expected,  the 
value  of  (A^-jV,)  for  DMDA-8320  is  much  smaller  than  others.  Its  value,  however,  is 
apparently  decreasing  with  the  shear  rate  unlike  all  others.  Consequently,  (^-NP;,)  of 
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DMDA-8320  show  stronger  shear  thinning  behavior  than  others.  The  reason  for  the 
crossover  between  HS-7028  and  DMDA-8320  is  unknown.  It  may  be  possible,  though, 
that  the  first  data  point  for  HS-7028  may  not  be  reliable  since  its  low  value  may  be  close 
to  the  low  limit  of  the  transducer  of  the  instrument. 

4.3  Coextrusion  Fiber  Spinning 

The  experimental  set-up  for  the  two-phase  fiber  is  shown  schematically  in  Figure 
24.  The  core  and  the  skin  materials  were  fed  by  two  19  mm  single  screw  extruders  into  a 
coextrusion  die  with  a 5.21  mm  diameter  hole.  The  schematic  of  the  coextrusion  die  is 
shown  in  Figure  25.  To  prevent  interfacial  instability  which  may  result  from  the  shear 
viscosity  mismatch  between  the  core  and  the  skin  materials,  the  feed  port  for  the  core 
material  was  made  to  extend  all  the  way  to  the  exit  of  the  coextrusion  die.  Consequently, 
the  core  and  the  skin  material  meet  right  at  the  exit  of  the  die  where  both  materials  are 
exposed  to  an  extensional  flow.  The  concentric  extrudate  was  then  pulled  under  tension 
and  quenched  as  it  passed  through  a water  bath  which  was  located  directly  below  the  die. 
The  melt  temperatures  of  the  core  and  skin  materials  were  maintained  at  200°C  and 
145°C,  respectively.  The  total  extrusion  output  was  held  constant  at  250  g/hr  throughout 
the  experiment  while  the  skin-layer  flow  rate  fraction  / was  varied  at  0.0,  0.2,  0.5,  0.8, 
and  1.0.  For  the  single-component  spinning  (i.e.,yN3.0  and  1.0),  the  same  material  (either 
LLDPE  or  LDPE)  was  fed  to  both  extruders  rather  than  using  only  one  extruder.  For 
each  experiment,  the  draw  span  was  also  varied  between  3 cm  and  10  cm.  In  the  draw 
span  the  polymer  melt  was  exposed  to  the  ambient  air  without  any  special  care  to 
maintain  the  melt  temperature.  The  spinning  process,  however,  was  presumably  at  a near 
isothermal  condition  since  the  radii  of  both  extrusion  die  and  the  fiber  were  very  large 
compared  to  those  of  existing  studies. 
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Core  material 


Figure  24.  Schematic  of  coextrusion  fiber  spinning  system. 
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Core  material 


—skin  material 
3.18mmID 


5.21mm 


skin  material 


Figure  25.  Schematic  of  coextrusion  die. 
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The  experiment  was  started  at  a low  take-up  speed  at  which  fibers  with  a uniform 
diameter  were  made.  The  take-up  speed  was  then  increased  gradually  in  a stepwise 
manner  while  taking  samples  at  each  speed.  A video  camera  which  was  connected  to  a 
TV  monitor  was  used  to  observe  the  variation  of  threadline  diameter.  In  all  cases  except 
the  single-component  LDPE  spinning,  the  draw  resonance  instability  became  visible 
beyond  a certain  critical  take-up  speed.  Most  experiments  were  continued  to  a much 
higher  take-up  speed  than  the  critical  value  as  long  as  fiber  breakage  did  not  occur.  The 
fiber  diameter  of  each  sample  was  then  measured  at  a 2.5  cm  interval  over  a 5 m distance 
to  determine  the  critical  draw  ratio,  the  maximum-to-minimum  diameter  ratio  as  well  as 
the  wave  length  of  the  instability.  Whenever  the  draw  resonance  was  visible,  its  period 
was  also  measured  using  a stop  watch  as  a supplementary  measurement. 

4.4  Results  and  Discussion 

The  variation  of  the  fiber  diameter  with  its  axial  distance  is  given  for  increasing 
values  of  draw  ratio  in  Figure  26.  This  figure,  which  is  for  f=0.23  (23%LDPE(Alathon- 
1540)/77%LLDPE(HS-7028))  at  a draw  span  of  3 cm  and  at  a flow  rate  of  250  g/hr, 
represents  the  general  trend  of  all  fiber  samples  excluding  the  one  for/=1.0  which  did  not 
show  draw  resonance.  It  is  apparent  from  the  figure  that  the  critical  draw  ratio  for  this 
experiment  is  between  51  and  57.  This  figure  also  indicates  that  the  draw  resonance  is 
sustained  beyond  the  critical  draw  ratio.  With  increasing  Dr,  both  amplitude  and  wave 
length  are  observed  to  be  increasing.  However,  the  period  of  the  instability  that  is 
determined  by  dividing  the  wavelength  by  the  take-up  speed  appears  to  remain  constant 
within  the  measurement  error.  A similar  behavior  was  observed  with  a polystyrene  by 
Matsumoto  and  Bogue  (1978).  In  the  experiment  of  Demay  and  Agassant  (1985)  which 
used  a linear  polyester,  however,  the  period  increased  slightly  with  draw  ratio. 


(xntn)Q  (raxn)d  (Tnra)a  (rara)d 
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Figure  26.  Variation  of  fiber  diameter  with  increasing  draw  ratio; 

(^=0.23,  draw  span  L= 3 cm). 
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In  Figure  27,  the  period  of  draw  resonance  instability  (/*)  is  plotted  as  a function 
of  L/w0  for  various  values  of / L and  w0  are  the  draw  span  and  the  average  extrusion 
velocity,  respectively.  It  seems  that  the  period  increases  linearly  with  L/w0  while  it 
decreases  as  the  thickness  of  the  more  viscoelastic  cladding  is  increased.  Also  included  in 
the  figure  is  the  prediction  of  the  linear  stability  theory  for  an  isothermal  Newtonian  fluid 
(solid  line  in  Figure  27)  in  which  the  imaginary  part  of  the  first  eigenvalue  is  14.02. 
From  Equation  (3.10),  t =0.448;  that  is 


/•  =0.448—  (4.2) 

w0 

It  appears  that  this  theoretical  prediction  is  in  good  agreement  with  the  single-phase 
LLDPE  data  for  which  the  linear  regression  yields 


0.445 — +1.5 


(4.3) 


This  linear  relation  between  t*  and  L/w0  is  presented  by  the  dashed  line  in  Figure  27.  The 
coefficient  of  correlation  for  the  regression  is  0.97.  In  calculating  the  average  extrusion 
velocity  w0,  we  have  used  the  cross-sectional  area  of  the  die  hole  as  the  basis.  If  the 
maximum  cross-sectional  area  of  the  fiber  at  the  die  swell  region  is  used  as  the  basis,  the 
corrected  w0  will  become  smaller  and  all  data  points  in  Figure  27  will  shift  to  the  right 
resulting  in  an  improved  agreement  between  the  LLDPE  data  and  the  Newtonian  theory. 
The  LLDPE,  which  shows  a shear  thinning  behavior,  is  by  no  means  a Newtonian  fluid. 
Nevertheless,  considering  the  fact  that  it  does  not  exhibit  a strain  hardening  behavior. 
The  agreement  between  the  experimental  data  and  the  Newtonian  theory  appears  to  be 
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Figure  27.  Period  f of  draw  resonance  versus  L/w0  for  the  coextrusion 

of  (a)  LLDPE(HS-7028)/LDPE(D YNH- 1 ):  (o)  100%  HS- 
7028,  ( — ) Newtonian  theory,  (□)  20%  DYNH-1,  (a)  50% 
DYNH-1,  (*)  80%DYNH-1;  (b)  LLDPE(HS-7028)/LDPE 
(Alathon-1540):  (o)100%HS-7028,  ( — )Newtonian  theory, 
(□)  20%  Alathon-1540,  (a)  50%  Alathon-1540,  (*)  80% 
Alathon-1540;  (c)  LLDPE(HS-7028)/LDPE:  (o)  100% 
HS-7028,  (— ) Newtonian  theory,  (□)  50%  DYNH-1,  (a) 
50%  Alathon-1540. 


Figure  27.  (continued) 


Figure  27.  (continued) 
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too  good  to  be  a simple  coincidence.  In  addition,  the  measured  critical  draw  ratio  was  18 
which  is  not  much  different  from  the  Newtonian  value  of  20.21.  The  critical  draw  ratio 
will  be  even  closer  to  20.21  if  the  maximum  cross-sectional  area  at  the  die  swell  region  is 
used  to  calculate  u0.  Furthermore,  in  our  experiment,  the  diameter  of  the  die  hole  (thus, 
the  diameter  of  the  fiber)  was  very  large  compared  with  those  of  other  experimental 
studies.  Consequently,  cooling  of  the  fiber  by  ambient  air  during  spinning  could  be 
relatively  small  creating  a near-isothermal  spinning  condition. 

The  decreasing  trend  of  period  with  increasing  skin-layer  flow  rate  fraction  / 
depicted  in  Figure  27  also  appears  to  agree  with  the  theoretical  prediction  described  in 
Figure  17.  In  the  present  experiment  the  critical  draw  ratio  increased  as  the  thickness  of 
the  skin-layer  material  was  increased  (Figure  30  and  31).  Therefore,  a large  value  of / 
corresponds  to  a large  critical  draw  ratio.  According  to  Figure  17,  the  draw  resonance 
occurring  at  a higher  draw  ratio  may  have  a smaller  period. 

Alathon- 1 540  among  the  two  different  skin  materials  is  expected  to  have  a larger 
relaxation  time  than  the  other.  According  to  the  theoretical  prediction,  then,  the  HS- 
7028/Alathon-1540  combination  should  have  a larger  critical  draw  ratio  and  a smaller 
period  when  the  instability  is  encountered  than  the  HS-7028/DYNH-1  combination.  The 
current  observations  seem  to  be  in  fact  in  accordance  with  the  predictions.  In  Figure  27c, 
the  periods  of  the  draw  resonance  instability  of  the  two  combinations  are  compared  when 
0.5.  The  HS-7028/Alathon-1540  combination  always  showed  smaller  period  than  the 
other.  Furthermore,  the  slope  of  the  apparently  straight  line  connecting  the  three  data 
points  is  also  smaller  than  the  other,  indicating  that  the  dimensionless  period  for  the  HS- 
7028/Alathon-1540  combination  is  smaller  than  the  other. 

In  Figure  26,  it  was  described  that  the  periodic  instability  was  sustained  beyond 
the  critical  draw  ratio  although  the  amplitude  (or  the  maximum-to-minimum  diameter 
ratio)  apparently  increased  with  the  draw  ratio.  The  maximum-to-minimum  diameter 
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ratio  of  the  unstable  fiber  is  plotted  as  a function  of  draw  ratio  for  20%-LDPE(DYNH-l) 
coextrusion  fiber  in  Figure  28.  It  is  apparent  from  the  figure  that  the  critical  draw  ratio 
for  this  experiment  is  in  the  vicinity  of  about  30.  Furthermore,  the  amplitude  of  the 
instability  (or  the  maximum-to-minimum  diameter  ratio)  appears  to  be  decreasing  with 
increasing  draw  span  although  the  critical  draw  ratio  remains  unchanged.  A similar  trend 
was  observed  with  a Newtonian  silicone  oil  by  Donnelly  and  Weinberger  (1975).  In  the 
experiment  of  Demay  and  Agassant  (1985),  however,  somewhat  different  trend  was 
observed  with  a linear  polyester  in  that  the  critical  draw  ratio  increased  as  the  draw  span 
increased  from  2 to  6 cm. 

The  variation  of  the  maximum-to-minimum  diameter  ratio  is  plotted  in  Figure 
29a  and  29b  for  three  different  values  of  / for  the  two  different  combinations, 
respectively.  The  changes  in  the  amplitude  of  the  instability  with  the  draw  span  for  all 
three  combinations  (f=Q. 2,  0.5  and  0.8)  also  show  the  same  trend  as  in  Figure  28.  When 
comparing  the  data  for  the  three  different  values  of / at  the  same  draw  span,  the  rate  of 
increase  of  the  maximum-to-minimum  diameter  ratio  with  the  draw  ratio  is  apparently 
slower  when  the  skin  flow  rate  ratio /is  increased.  This  may  be  another  manifestation  of 
the  stabilizing  effect  of  the  skin  layer.  These  figures  also  imply  that  the  critical  draw  ratio 
increases  with  the  skin  flow  rate  fraction / This  behavior  is  summarized  in  Figure  30. 

In  Figure  30,  the  critical  draw  ratio  is  plotted  as  a function  of  the  skin-layer  flow 
rate  fraction/for  the  two  different  material  combinations.  The  mode  of  instability  for  all 
data  was  draw  resonance  except  for  the  single-phase  spinning  of  LDPE  (i.e.,  ./HO). 
When^l.0,  draw  resonance  was  never  observed.  Instead,  the  fiber  snapped  at  a relatively 
large  draw  ratio.  It  should  be  noted  that  the  "maximum"  draw  ratio  before  the  fiber 
breakage  increased  with  increasing  draw  span  whereas  the  "critical"  draw  ratio  for  all 
other  fiber  compositions  was  insensitive  to  the  draw  span.  A similar  behavior  was  also 
observed  in  the  single-phase  experiment  of  Chang  and  Denn  (1979). 


niax/ min  diameter  ratio 
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Figure  28.  The  maximum-to-minimum  diameter  ratio  as  a function  of 
draw  ratio  for  the  combination  of  20%-LDPE(DYNH- 
1 )/80%-LLDPE(HS-7028)  at  three  different  draw  spans:  ( 
o)  L= 3 cm;  (□)  L= 6 cm;  (a)  L=  10  cm. 


max/ min  diameter  ratio 
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Draw  ratio 


Figure  29.  Variation  of  the  maximum-to-minimum  diameter  ratio  for 

three  combinations  of  (a)  LDPE(DYNH-1)/LLDPE(HS- 
7028);  (b)  LDPE(Alathon-1540)/LLDPE(HS-7028)  at 
different  draw  spans:  (o)  L= 3 cm;  (□)  L= 6 cm;  (a)  L=10 
cm. 


max/ min  diameter  ratio 


94 


Draw  ratio 


Figure  29.  (continued) 
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Figure  30.  Critical  draw  ratio  as  a function  of  skin-layer  flow  rate 

fraction  for  different  combinations  at  different  draw  spans: 
(o)  L= 3 cm;  (□)  1=6  cm;  (a)  1=10  cm.  Core:  LLDPE(HS- 
7028);  Skin:  LDPE  ( — : Alathon-1540, : DYNH-1). 
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To  this  point,  the  draw  resonance  characteristics  of  only  two  different  material 
combinations  have  been  discussed  in  which  the  core  material  was  fixed  with  HS-7028. 
Similar  experiments  have  been  also  conducted  with  the  DMDA-8320  as  the  core  whose 
shear  viscosity  is  much  smaller  than  that  of  HS-7028.  Thus,  the  following  four  different 
combinations  were  investigated  in  this  experimental  study: 

( 1 ) LLDPE(HS-7028)/LDPE( Alathon- 1 540) 

(2)  LLDPE(DMDA-8320)/LDPE( Alathon- 1 540) 

(3)  LLDPE(HS-7028)/LDPE(DYNH- 1 ) 

(4)  LLDPE(DMDA-8320)/LDPE(DYNH-l) 

The  results  for  (1)  and  (2)  have  been  already  discussed  extensively.  The  difference 
between  (1)  and  (2)  (or  between  (3)  and  (4))  are  the  shear  viscosity  of  the  core  material. 
When  compared  with  the  theoretical  analysis,  the  combination  (1)  has  a higher  value  of 
M (=fic  / ns)  than  combiantion  (2).  Equivalently,  (3)  has  a higher  M than  (4). 

According  to  the  theoretical  predicition  described  in  Figure  19,  when  M is 
decreased,  the  relative  influence  of  the  viscoelastic  skin  layer  increases  resulting  in  a 
slight  increase  in  critical  draw  ratio  Dr*  if  draw  resonance  is  still  encountered  or  a large 
increase  in  Dr*  if  the  Deborah  number  happens  to  be  on  the  right  hand  side  of  the 
unstable  nose-shaped  envelope.  Unfortunately  neither  of  these  were  observed  in  the 
present  experiment.  The  measured  dependence  of  the  critical  draw  ratio  on  the  skin-layer 
flow  rate  fraction /Hs  given  in  Figure  31  for  all  four  combinations.  In  all  coextrusion 
experiments,  the  mode  of  instability  was  draw  resonance.  When  (1)  and  (3)  are 
compared,  stronger  influence  of  Alathon- 1540  than  DYNH-1  is  apparent.  When  (2)  and 
(4)  are  compared,  however,  no  apparent  difference  between  Alathon- 1540  and  DYNH-1 
is  observed.  This  result  is  rather  perplexing  since  the  viscoelasticity  effect  of  the  skin 
layer  should  be  more  pronounced  when  the  shear  viscosity  of  the  core  is  smaller. 
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Figure  31.  Critical  draw  ratio  as  a function  of  skin-layer  flow  rate 

fraction  for  different  combinations  at  L= 6 cm: 

(1)  Core:LLDPE(HS-7028);  Skin:LDPE(Alathon-1540) 

(2)  Core:LLDPE(DMDA-8320);Skin:LDPE(Alathon- 1 540) 

(3)  Core:LLDPE(HS-7028);  Skin:LDPE(DYNH-l) 

(4)  Core:LLDPE(DMDA-8320);Skin:LDPE(DYNH-l) 
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The  theoretical  prediction  that  the  highest  draw  ratio  for  a stable  coextrusion 
spinning  can  be  larger  than  that  for  the  single-phase  spinning  of  either  fluid  was  not  also 
observed  in  the  present  experiments.  However,  it  does  not  necessarily  mean  that  the 
theory  is  incorrect.  In  all  coextrusion  experiments  with  the  four  different  combination  of 
materials,  the  mode  of  instability  was  draw  resonance.  This  may  indicate  that  the 

viscoelasticity  of  the  skin  materials  are  not  large  enough  to  exhibit  the  theoretical 
prediction. 


4.5  Summary  and  Conclusion 

Draw  resonance  instability  in  bicomponent  fiber  spinning  has  been  studied 
experimentally  using  four  polymers  with  different  rheological  properties;  two  linear  low- 
density  polyethylenes  (LLDPE)  and  two  branched  low-density  polyethylenes  (LDPE). 
The  LLDPEs  which  are  very  prone  to  draw  resonance  were  used  as  the  core  material, 
whereas  the  LDPEs  which  do  not  show  draw  resonance  instability  were  used  as  the 
cladding  material.  The  results  indicate  that 

1)  In  all  coextrusion  experiments  with  the  given  materials,  the  draw  resonance  was 
always  the  dominant  mode  of  instability  even  when  the  LDPE  fraction  was  as  high 
as  0.8.  The  critical  draw  ratio,  however,  increased  substantially  with  increasing 
LDPE  fraction. 

2)  Draw  resonance  instability  was  sustained  beyond  the  critical  draw  ratio.  The 
amplitude  of  the  instability  was  observed  to  be  increasing  with  the  draw  ratio 
although  its  period  remained  unchanged. 

3)  The  critical  draw  ratio  appeared  to  be  independent  of  draw  span  although  the 
maximum-to-minimum  diameter  ratio  decreased  with  increasing  draw  span. 

4)  The  period  of  the  instability  increased  linearly  with  L/w0  while  it  decreased  with 
increasing  LDPE  content. 
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Many  of  the  experimental  results  are  in  reasonable  agreement  with  theoretical 
predictions.  However,  one  of  the  most  interesting  predictions  regarding  a kind  of 
synergism  in  terms  of  maximum  draw  ratio  was  not  observed  with  the  given  choice  of 
materials.  Although  the  present  results  are  specifically  for  a fiber  spinning  process,  the 

general  trend  should  be  also  relevant  to  other  fast-drawing  processes  such  as  extrusion 
coating  or  film  casting. 


CHAPTER  5 

FABRICATION  OF  PLASTIC  OPTICAL  FIBERS 


5.1  Introduction 

A direct  application  of  the  two-phase  fiber  spinning  process  may  be  in  the 
manufacture  of  plastic  optical  fibers.  In  this  Chapter,  a tiew  processing  method  for  the 
fabrication  of  plastic  optical  fibers  is  investigated  which  is  a modified  coextrusion 
process.  The  new  method  has  been  developed  based  on  the  fundamental  understanding  of 
the  two-phase  fiber  spinning  process  which  has  been  investigated  in  the  previous 
Chapters. 

Plastic  optical  fibers  offer  a number  of  advantages  compared  to  silica  glass 
optical  fibers.  Besides  sizable  cost  savings  in  materials,  toughness  and  durability  of 
plastics  allow  plastic  optical  fibers  to  be  handled  without  special  care.  In  addition,  the 
mechanical  flexibility  of  plastic  allows  these  fibers  to  have  large  size  cores  that  may 
range  from  0.5  to  1.0  mm  in  diameter.  Consequently,  a low-cost,  large-area  light  source 
can  be  used  and  the  connection  between  the  fibers  becomes  easier  and  less  costly. 
Furthermore,  a high  refractive-index  difference  can  be  achieved  between  core  and 
cladding  materials  yielding  a higher  numerical  aperture  than  glass  optical  fibers. 

Plastic  optical  fibers,  however,  have  much  larger  attenuation  of  light  signals  than 
glass  fibers  due  to  their  large  intrinsic  absorption  loss.  In  case  of  PMMA 
(poly(methylmethacrylate))  optical  fibers,  the  theoretical  lower  loss  limit  is  about  105 
dB/km  at  650  nm  and  the  commercial  PMMA  fibers  have  a typical  attenuation  of  about 
200-400  dB/km  at  the  same  wave  length.  Silica  glass  optical  fibers,  on  the  other  hand, 
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are  readily  manufactured  these  days  that  have  an  attenuation  as  small  as  0.5  dB/km  at 
1.30pm.  Despite  the  high  attenuation,  plastic  optical  fibers  are  economically  very 

attractive  and  are  better  suited  than  glass  fibers  to  many  applications  such  as  fiber  optic 
sensors,  local  area  network,  automotive  and  medical  applications. 

Commercial  plastic  optical  fibers  are  typically  made  of  a polystyrene(PS)  core 
and  a PMMA  cladding  or  a PMMA  core  and  a fluoropolymer  cladding.  These  materials 
are  amorphous  thermoplastics  that  are  optically  transparent  and  drawable  to  fibers.  In  the 
fabrication  of  these  plastic  optical  fibers,  a cylindrical  rod  (or  a preform)  is  made  first 
using  highly  purified  monomers.  The  preform  is  then  heated  and  drawn  to  optical  fibers. 
This  two-step  process  which  may  be  called  as  "thermal-drawing"  process  is  similar  to  the 
vapor  phase  oxidation  process  for  glass  optical  fibers  as  they  both  make  preforms  and 
then  heat-draw  fibers. 

Since  PS,  PMMA,  and  fluoropolymer  cladding  are  thermoplastic,  a coextrusion 
process  may  be  also  applicable  as  an  alternative  method  to  fabricate  optical  fibers.  In  a 
coextrusion  process  the  core  and  the  cladding  materials  are  extruded  simultaneously  into 
a coextrusion  die  using  two  separate  extruders,  and  a clad  fiber  (i.e.,  optical  fiber)  is 
drawn  out  of  the  coextrusion  die  to  a desired  diameter.  This  process  is,  in  a sense, 
similar  to  the  direct-melt  method  (or  double-crucible  method)  for  glass  fibers  in  which 
optical  fibers  are  made  directly  from  the  molten  state  of  silicate  glasses  following  a 
traditional  glass-making  procedure. 

The  coextrusion  method  is  a single-step  process  that  may  appear  to  be  less 
complicated  than  the  heat-drawing  two-step  process.  This  method,  however,  is 
commercially  less  prevalent  due  to  its  susceptibility  to  contamination  by  impurities.  In 
the  coextrusion  process,  polymers  in  a pellet  form  are  used  for  both  core  and  cladding. 
Although  extra  care  may  be  taken  to  maintain  the  high  purity  of  the  materials,  the  pellets 
are  highly  susceptible  to  contamination  by  foreign  materials  such  as  dust  particles  during 
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the  materials  handling  (e.g.,  transportation  and  feeding  for  extrusion).  In  the  two-step 
thermal-drawing  process,  on  the  other  hand,  a purified  core  material  may  never  be 
exposed  to  foreign  substances  throughout  the  process  from  polymerization  to  heat- 
drawing resulting  in  superior  quality  optical  fibers. 

In  these  conventional  methods  (thermal-drawing  or  coextrusion),  the  materials 
that  can  be  used  to  make  optical  fibers  are  limited  to  thermoplastics.  Non-thermoplastic 
materials  cannot  be  fabricated  into  optical  fibers  despite  their  excellent  optical  properties. 
To  fabricate  optical  fibers  using  various  classes  of  materials,  a new  processing  method 
has  been  developed.  The  new  method  is  a modified  coextrusion  process.  Unlike  the 
conventional  coextrusion,  however,  only  the  cladding  material  has  to  be  thermoplastic  in 
this  process  and  the  core  material  can  be  any  class  of  materials;  liquids  (monomers  or 
prepolymers),  elastomeric  gels,  thermoplastic  or  thermosetting  materials. 

The  process  modeling  study  described  in  the  previous  Chapters  suggests  that  the 
mechanics  of  the  two-phase  flow  is  dictated  by  the  more  viscoelastic  fluid  of  the  two 
components  under  typical  processing  conditions.  Thus,  in  the  new  method,  several 
different  types  of  materials  can  be  used  as  the  core  as  long  as  the  skin  material  has  high 
enough  viscoelasticity  to  control  the  overall  flow  dynamics.  In  the  following  sections,  a 
brief  description  of  the  new  method  and  the  experimental  results  on  the  feasibility  of  the 
new  method  are  discussed  including  the  advantages  and  limitations  of  the  new  method. 

5.2  New  Processing 

In  Figure  32,  a schematic  of  the  new  processing  method  is  given  in  which  a 
thermoplastic  cladding  material  is  spun  into  a hollow  fiber  while  a core  material  is  fed  by 
a piston-cylinder  unit  into  the  core  of  the  hollow  fiber  at  the  same  time.  In  this  process, 
the  cladding  material  has  to  be  a thermoplastic  as  it  is  processed  by  an  extruder.  The 
core  material,  however,  can  be  any  type  of  materials  which  may  be  fed  by  an  extruder,  a 


Figure  32.  Experimental  set-up  for  the  spinning  of  plastic  optical  fibers 
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piston-cylinder  assembly,  or  a gear  pump  depending  upon  its  rheological  properties. 
After  fibers  are  made  by  this  process,  the  core  material  may  be  either  cross-linked  or 
further  polymerized  (if  the  core  is  a low  molecular  weight  prepolymer)  to  enhance  the 
mechanical  strength  of  the  fiber  if  such  post  processing  is  necessary.  In  this  way,  the 
process  of  making  preforms  can  be  eliminated  and  both  thermosetting  and  thermoplastic 
materials  can  be  used  for  optical  fibers.  In  addition,  contamination  by  impurities 
(including  dust)  can  be  minimized  since  purified  liquids  (monomers  or  prepolymers)  are 
fed  directly  into  the  core  without  contacting  any  foreign  substances  or  without  being 
exposed  to  ambient  air. 

While  the  most  important  feature  of  the  new  method  is  its  applicability  to  various 
materials  including  those  which  cannot  be  made  into  optical  fibers  by  the  conventional 
methods,  the  new  method  also  has  some  limitations.  First,  the  liquid  phase  core  should 
be  a non-solvent  of  the  thermoplastic  cladding.  Otherwise,  the  cladding  will  be  seriously 
damaged  before  the  core  material  is  completely  reacted.  Furthermore,  the  polymerization 
or  cross-linking  reaction  should  take  place  at  a temperature  below  the  glass  transition 
temperature  of  the  cladding  material  (or  below  the  melt  temperature  if  the  cladding  is  a 
semi-crystalline  polymer).  Finally,  the  density  variation  involved  in  the  phase  change  of 
the  core  material  should  be  minimal.  Otherwise,  the  fibers  may  deform  or  bubbles  may 
be  formed  during  the  post-reaction.  Thus,  prepolymers  with  a certain  level  of  reaction 
are  preferable  to  monomers  as  a core  material.  These  appear  to  be  rather  severe 
restriction.  Nevertheless,  many  different  combinations  of  materials  are  possible  with 
which  the  restrictive  conditions  can  be  met. 
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5.3  Optical  Fiber  Spinning 

5.3.1  Experiments 

A schematic  of  the  equipmental  set-up  is  given  in  Figure  32.  A thermoplastic 
cladding  material  is  fed  into  a coextrusion  die  by  a 3/4"  extruder  while  a core  material  is 
fed  by  a piston-cylinder  assembly  which  has  been  custom-made.  The  volume  of  the 
cylinder  is  about  100  cm3  which  can  effectively  handle  small  quantity  of  materials  to 
perform  pilot  scale  experiments.  The  lead  screw  which  drives  the  piston  is  driven  by  a 
motor  which  provides  an  accurate  feed  volume  control  of  core  materials  from  4 cm3/min 
to  40  cm3/min. 

The  design  of  the  coextrusion  die  is  given  schematically  in  Figure  25.  The  most 
important  features  of  this  design  are  (1)  the  protrusion  of  the  core  capillary  to  the  end  of 
the  die  and  (2)  a short  thermal  contact  time  between  core  and  cladding  materials.  To 
handle  a large  viscosity  difference  between  the  core  and  cladding  materials,  the  capillary 
that  feeds  the  core  material  needs  to  be  stretched  to  the  end  of  the  die  exit.  When  a 
crosslinkable  core  material  is  used,  excessive  heating  should  be  avoided  to  repress  the 
reaction  in  the  coextrusion  die.  In  the  current  die  design,  the  capillary  diameter  is  small 
(3.21  mm)  and  its  length  short  (7.5  cm).  Thus,  the  residence  time  of  the  core  material  is 
relatively  small  minimizing  its  thermal  contact  with  the  hot  environment  which  may  be 
kept  at  150  to  170°C  depending  upon  the  choice  of  the  cladding  material. 

Using  the  prescribed  experimental  set-up,  spinning  experiments  have  been 
conducted  with  two  different  types  of  core  materials;  (1)  styrene-hexyl  methacrylate 
copolymer  and  (2)  polysiloxanes.  The  most  important  properties  of  the  optical  polymers 
are  their  transparency  and  refractive  indices  to  achieve  a minimum  level  of  light  signal 
loss.  Since  optical  polymers  are  highly  proprietary,  the  raw  materials  are  very  difficult  to 
obtain.  Consequently,  the  choice  of  materials  was  limited  in  the  present  experiment.  In 
any  case,  the  major  emphasis  of  the  experiments  was  on  the  processing  aspect  of  the  new 
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method  to  evaluate  its  feasibility  to  fabricate  physically  sound  optical  fibers  that  are  free 
from  critical  defects  (e.g.,  bubbles  or  diameter  variation),  and  relatively  little  attention 
has  been  given  to  measure  the  attenuation  length  of  the  fiber.  To  achieve  a high 
attenuation  length,  extensive  purification  of  the  core  materials  should  be  preceded  which 
is  not  necessary  for  the  present  objective. 

5.3.2  Styrene  Core  Fibers 

Styrene  is  one  of  the  most  common  materials  that  are  used  as  the  core  of  plastic 
optical  fibers,  and  it  has  been  chosen  as  the  first  material  for  the  feasibility  study  of  the 
new  processing  method.  Purified  styrene  monomer  was  prepolymerized  in  a flask  at 
120°C  for  three  hours.  The  prepolymer  was  then  poured  into  the  heated  cylinder  at 
120°C  where  it  was  then  kept  for  two  more  hours  for  further  reaction  (Figure  32).  Since 
the  density  increase  during  the  polymerization  is  rather  substantial  (about  1 6%),  a high 

conversion  of  polymerization  reaction  was  necessary  before  the  coextrusion  spinning  to 
avoid  bubble  formation  during  the  post-polymerization.  The  core  material  was  then  fed 
into  the  coextrusion  die  by  a piston  to  fabricate  1 mm  optical  fibers.  For  this  experiment 
an  ethylene-vinyl  acetate  copolymer  was  used  as  a cladding  material.  The  conventional 
PMMA  was  not  used  since  the  core  material  which  might  still  contain  styrene  monomer 
could  seriously  damage  the  PMMA  cladding.  After  the  spinning,  the  fiber  was  kept  in  an 
oven  at  100°C  for  24  hours  to  complete  the  polymerization  of  the  core  material.  Due  to 

the  high  refractive  index  of  the  cladding  material  (n=1.50),  the  fiber  was  not  expected  to 
have  a high  attenuation  length.  Nevertheless,  this  experiment  demonstrated  the 
feasibility  of  the  new  processing  method  in  which  a polymerization  reaction  can  be 
involved. 
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5.3.3  Modification  of  the  Spinning  Equipment 

The  piston-cylinder  unit  in  Figure  32  served  both  as  a reactor  and  as  a feeder. 
One  limitation  of  this  unit  as  a reactor  was  its  inability  to  provide  good  mixing  during 
polymerization  reaction  resulting  in  a poor  molecular  weight  control.  Consequently,  the 
polystyrene  core  of  the  fiber  was  very  brittle.  Good  mixing  is  especially  important  for 
the  copolymerization  which  is  intended  to  fabricate  flexible  optical  fibers  by 
incorporating  n-hexyl  methacrylate  into  styrene.  Thus,  the  experimental  set-up  has  been 
modified  as  described  in  Figure  33.  With  this  modification,  not  only  the  control  of  the 
polymerization  reaction  was  more  efficient  with  uniform  mixing,  but  also  the  spinning 
experiment  could  be  conducted  for  a longer  period  time  due  to  the  increase  in  reactor 
volume  ( 1 000  vs.  1 00  ml). 

5.3.4  Stvrene-Hexvl  Methacrylate  Copolymer 

Spinning  experiments  were  conducted  using  the  modified  experimental  set-up 
with  styrene  and  n-hexyl  methacrylate  copolymer  as  the  core  material.  Copolymerization 
of  styrene  and  n-hexyl  methacrylate  is  to  obtain  a flexible  material  whose  glass  transition 
temperature  is  below  the  room  temperature  so  that  a large  diameter  (greater  than  1mm) 
plastic  optical  fibers  can  be  manufactured.  The  molar  ratio  of  styrene  and  n-hexyl 
methacrylate  was  varied  at  80/20  and  50/50,  respectively.  The  refractive  indices  of 
polystyrene  and  poly(n-hexyl  methacrylate)  are  1.59  and  1.48,  respectively  and  the 
refractive  index  of  the  copolymer  is  expected  to  be  linearly  dependent  on  the  mole 
fraction  of  n-hexyl  methacrylate  between  the  two  values.  Therefore,  as  the  molar  content 
of  n-hexyl  methacrylate  is  increased,  the  refractive  index  of  the  core  decreases  whereas 
the  fiber  becomes  more  flexible. 

In  these  experiments  n-butyl  mercaptan  was  added  at  20  mmol//  as  a chain 
transfer  agent  and  the  polymerization  temperature  was  maintained  at  120°C.  Fiber 
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spinning  was  initiated  after  four  hours  of  reaction  using  the  same  ethylene-vinyl  acetate 
copolymer  as  the  cladding  material.  Fibers  with  1. 5-2.0  mm  diameter  were  fabricated 
and  aged  in  an  oven  at  100°C  for  post-polymerization.  After  the  aging  in  the  oven,  it 
was  found  that  the  cross-section  of  the  fiber  was  rather  elliptic  than  circular.  This 
deformation  is  presumably  due  to  the  volume  change  resulting  from  the  loss  of 
monomers  by  diffusion  during  the  aging.  This  result  suggests  that  further  polymerization 
is  necessary  before  the  spinning  experiment.  Further  polymerization,  however,  is  not 
desirable  with  the  current  reactor  set-up  since  the  viscosity  of  the  prepolymer  will  be 
increased  substantially  causing  cavitation  problem  during  pumping.  This  problem, 
however,  can  be  resolve  if  the  current  reactor  is  replaced  by  a pressurizable  one  so  that 
the  pressure  can  be  raised  at  the  feeding  port  of  the  gear  pump. 

5.3.5  Polvsiloxanes 

Polysiloxanes  based  on  the  monomers  of  dimethyl,  methyl-phenyl,  and  diphenyl 
siloxane  are  optically  transparent.  It  is  known  that  the  refractive  index  (n)  of 
polysiloxane  increases  linearly  with  its  phenyl  content.  Among  various  polysiloxanes, 
poly(dimethyl  diphenylsiloxane)  and  poly(methylphenylsiloxane)  were  initially  used  as 
the  core  along  with  a PMMA  cladding.  The  poly(dimethyl  diphenylsiloxane)  was  an 
elastomeric  gel,  whereas  the  poly(methylphenylsiloxane)  was  a viscous  liquid.  Since  the 
viscosity  of  the  latter  was  so  low,  it  was  fed  by  a liquid  gear  pump  instead  of  the  piston- 
cylinder  unit.  For  both  cases  physically  intact  optical  fibers  could  be  made.  Since  both 
of  these  materials  were  chemically  stable,  crosslinking  reaction  could  not  be  induced  and 
the  core  materials  remained  as  a gel  and  as  a liquid  phase,  respectively. 

Among  many  different  commercial  polysiloxanes,  poly(dimethyl  siloxane)  with 
vinyl-phenyl-methyl  termination  was  the  only  crosslinkable  one  we  could  Find.  Although 
its  refractive  index  is  as  low  as  1.39  due  to  its  low  phenyl  content,  this  material  was  used 
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to  conduct  a spinning  experiment  which  involved  a crosslinking  reaction.  The  reaction 
study  in  a test  tube  has  indicated  that  the  volume  change  (or  density  change)  is  negligible 
during  the  crosslinking  reaction.  Thus,  unlike  the  styrene-based  material,  polysiloxane 
could  be  coextruded  as  a liquid  phase  core  followed  by  post  polymerization  after 
spinning.  For  the  spinning  experiment  with  this  core  material,  benzoyl  peroxide  was 
added  at  0.2%  as  an  initiator  of  the  crosslinking  reaction  and  an  ethylene-vinyl  acetate 
copolymer  was  used  as  the  cladding.  After  the  spinning  process  the  fibers  were  put  in  an 
oven  at  100°C  to  induce  the  crosslinking  reaction.  Unlike  the  reaction  study  in  test  tubes, 

the  crosslinking  reaction  was  found  to  be  rather  irregular  in  that  the  reaction  was 
complete  with  some  samples  while  it  was  not  with  others.  It  is  possible  that  the 
irregularity  may  result  from  a low  reaction  temperature,  low  initiator  content,  or  some 
unknown  interference  of  the  cladding  material.  Further  study  should  follow  to  resolve 
this  problem. 

The  polysiloxanes  were  specifically  studied  in  conjunction  with  the  US 
Department  of  Energy  Project  to  develop  "radiation  hard  scintillating  plastic  optical 
fibers"  for  the  Superconducting  Super  Collider  (SSC).  In  the  harsh  environment  of  the 
SSC,  the  commercial  optical  fibers  based  on  PS-PMMA  materials  are  inadequate  due  to 
severe  radiation  damage.  Siloxane  based  materials  were  identified  as  one  of  the  possible 
choices  for  the  radiation  hard  optical  fibers,  and  attempts  were  made  in  the  present  study 
to  convert  the  polysiloxane  into  optical  optical  fibers.  The  materials  those  have  been 
evaluated  are  either  in  liquid  or  gel  phase  and  consequently,  the  conventional  thermal- 
drawing  method  is  not  applicable.  In  this  study,  it  has  been  demonstrated  that  such 
materials  can  in  fact  be  fabricated  into  optical  fibers  by  the  new  processing  method. 
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5.3.6  Uniformity  of  Fiber  Radius 

Uniformity  of  fiber  radius  is  one  of  the  most  important  physical  aspects  which 
has  significant  influence  on  the  optical  property  of  the  optical  fibers.  Throughout  the 
experiments,  the  fiber  radius  was  not  only  monitored  visually  but  also  inspected 
physically  using  a micrometer  within  an  accuracy  of  0.01mm.  For  at  least  one-meter 
long  fibers  the  outermost  diameter  was  measured  at  a 10  cm  interval.  For  an  average 
diameter  of  1 mm,  the  standard  deviation  was  within  3%.  This  level  of  uniformity  may 
be  tolerable  for  the  purpose  of  optical  fibers. 

5.4  Summary  and  Conclusion 

A new  processing  method  has  been  developed  for  the  fabrication  of  plastic  optical  fibers 
in  which  a thermoplastic  cladding  material  is  spun  into  a hollow  fiber  while  a core 
material  is  fed  into  the  core  of  the  hollow  fiber  at  the  same  time.  In  conventional 
methods  both  core  and  cladding  materials  should  be  thermoplastics.  In  the  new  method, 
however,  only  the  cladding  material  has  to  be  a thermoplastic  and  the  core  material  can 
be  any  class  of  materials  including  thermosetting  polymers.  The  experimental  study  has 
focused  on  assessing  the  feasibility  of  the  new  fabrication  method  to  fabricate 
physically  sound  plastic  optical  fibers.  Lab  scale  experiments  have  demonstrated  that 
the  new  method  is  capable  of  fabricating  defect-free  clad  fibers  using  several  different 
types  of  materials.  The  new  method  is  possible  because  the  mechanics  of  the 
bicomponent  spinning  flow  is  controlled  by  the  viscoelastic  cladding  material  under 
typical  spinning  conditions  as  predicted  by  the  process  modeling  study. 


CHAPTER  6 

SUMMARY  AND  CONCLUSION 


The  two-phase  flow  of  a bicomponent  coextrusion  fiber  spinning  has  been  studied 
theoretically  and  experimentally.  A simple  model  has  been  adopted  for  the  theoretical 
analysis  in  which  a Newtonian  and  an  upper-convected  Maxwell  fluid  are  the  core  and 
the  skin  layers,  respectively.  This  model  was  chosen  to  investigate  a flow  situation  in 
which  the  two  fluids  have  quite  different  extensional  rheology.  The  results  indicate  that 
the  viscoelasticity  of  the  skin  layer  plays  a significant  role  in  determining  the  overall 
mechanics  of  the  two-phase  flow.  Thus,  the  stability  of  the  flow  is  also  significantly 
influenced  by  the  viscoelastic  skin  layer  and  thereby  the  onset  of  draw  resonance  of  a 
Newtonian  fluid  is  delayed  to  a higher  value  of  draw  ratio  in  the  two-phase  coextrusion 
flow.  This  stabilizing  effect  becomes  more  prominent  as  the  viscoelasticity  represented 
by  the  Deborah  number  is  increased.  As  a consequence  of  the  competing  influence  of  the 
purely  viscous  and  viscoelastic  fluid,  an  interesting  behavior  has  been  also  predicted  in 
that  the  highest  draw  ratio  for  a stable  coextrusion  spinning  can  be  larger  than  that  of  the 
single-phase  spinning  of  either  fluid  (Newtonian  or  UCM)  under  a certain  condition. 

In  the  experimental  study  linear  low-density  polyethylenes  (LLDPE)  and  low- 
density  polyethylenes  (LDPE)  were  used  as  the  core-  and  the  skin-layer  materials, 
respectively.  This  choice  of  materials  has  been  made  to  simulate  the  flow  situation 
considered  in  the  theoretical  investigation  as  closely  as  possible.  The  two  materials  are 
chemically  identical  although  their  molecular  structures  are  different.  Thus,  they  have 
good  chemical  compatibility  ensuring  the  continuity  of  velocity  at  the  interface. 
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However,  they  are  known  to  have  distinctly  different  extensional  rheology  resulting  from 
the  differences  in  their  molecular  structure.  While  the  LLDPEs  are  very  prone  to  draw 
resonance  instability,  the  LDPEs  are  not.  It  is  due  to  the  larger  relaxation  time  of  the 
LDPEs  compared  to  the  LLDPEs.  The  experimental  results  may  be  summarized  as 
follows: 

1)  The  critical  draw  ratio  for  the  onset  of  draw  resonance  increases  as  the  flow  rate 
fraction  of  the  viscoelastic  skin  layer  is  increased. 

2)  Draw  resonance  instability  appears  to  be  sustained  beyond  the  critical  draw  ratio. 
With  increasing  draw  ratio,  the  amplitude  of  the  instability  increases  whereas  its 
period  remains  unchanged. 

3)  The  critical  draw  ratio  is  apparently  independent  of  the  draw  span  although  the 
maximum-to-minimum  diameter  ratio  is  observed  to  be  decreasing  with  increasing 
draw  span. 

4)  The  period  of  the  instability  increases  linearly  with  L/wn,  the  ratio  of  the  draw  span 
to  the  average  velocity  at  the  die  exit. 

The  linear  stability  theory  provides  only  a limited  information  about  the 
dynamical  behavior  of  the  instability  for  small  time.  Nevertheless,  many  aspects  of  the 
experimental  observations  seem  to  be  in  accordance  with  the  predictions  of  the  linear 
stability  analysis  at  least  qualitatively.  The  prediction  that  the  highest  draw  ratio  for  a 
stable  coextrusion  spinning  can  be  larger  than  that  of  the  single-phase  spinning  of  either 
fluid,  however,  has  not  been  observed  with  the  current  choice  of  materials. 

As  a practical  application  of  the  two-phase  fiber  spinning  process,  a new 
processing  method  has  been  developed  to  fabricate  plastic  optical  fibers  using  various 
type  of  materials.  The  new  process  method  is  possible  since  the  more  viscoelastic 
material  of  the  two  components,  which  is  typically  the  skin-layer  material,  dictates  the 
overall  mechanics  of  the  flow  as  predicted  by  the  theoretical  model.  While  the 
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conventional  coextrusion  process  is  for  extruding  only  thermoplastic  materials,  current 
study  suggests  a way  to  explore  a variety  of  core  material  options  for  plastic  optical 
fibers. 

The  study  presented  in  this  thesis  is  specifically  for  a fiber  spinning  process.  The 
mathematical  similarities  in  the  governing  equations,  however,  suggest  that  the 
mathematical  approaches  and  the  results  regarding  the  characteristics  of  the  two-phase 
flow  are  also  relevant  to  other  fast-drawing  process  such  as  extrusion  coating  or  film 
casting.  The  present  model  is  very  simple  in  that  simplest  constitutive  models  have  been 
used  for  an  isothermal  flow  situation  in  the  absence  of  the  gravitational  force  or  the 
surface  forces.  Nevertheless,  it  is  believed  that  the  present  model  represents  the 
fundamental  and  the  essential  features  of  the  two-phase  flow  system  at  least  qualitatively, 
so  that  it  may  serve  as  a basis  for  more  rigorous  analyses  using  more  complicated  but 
more  realistic  constitutive  models. 


APPENDIX  A 

THE  FINITE  DIFFERENCE  APPROACH 
OF  EIGENVALUE  PROBLEMS 


Several  mathematical  approaches  can  be  taken  to  solve  the  eigenvalue  problem. 
In  the  present  analysis,  the  eigenvalue  As  are  determined  by  reducing  the  eigenvalue 

problem  to  an  algebraic  eigenvalue  problem  by  the  use  of  finite  differences. 

The  draw-down  region  was  discretized  into  n equal  spaces  and  the  three  point 
central  finite  difference  was  used  as  follows  for  the  first  and  the  second  order  derivatives: 


= y'+\  ~ X-i 
2 h 


(A.l) 
(A. 2) 


Here  y'  represents  all  the  variables  involved  in  the  fiber  spinning  flow  problem;  w\  a , 
T__  and  T\  h is  the  equal-step  size  which  is  1/n.  Using  equation  (A.l)  for  the  continuity 
equation  (3.5)  is  given  as 


ai+ 1 ai-\ 


2 h 


) = -M*; 


(A. 3) 


Here  i=l  and  n+1  correspond  to  the  two  boundary  points  at  z= 0 and  1,  respectively.  At 
the  initial  point  z=0  (i.e.,  i=l),  all  variables  are  specified  as 
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<=w;  = t:,=7;'=w^=o 


(A.4) 


when  i=n+l,  since  the  central  differencing  cannot  be  used,  the  backward  differencing 
was  used  as  follows: 
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(A. 5) 


Using  these  finite  difference  equations  for  all  nodal  points  within  the  boundary,  the 
continuity  equation  is  formulated  in  a matrix  form  as 
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The  velocity  at  z=l  is  fixed.  Thus,  H’*+1=0  and  consequently,  the  dimension  of  the 
unknow  vector  w*  is  (n-1)  whereas  that  of  a ’ is  n. 

Similary,  the  axial  force  balance  (3.6)  can  be  expressed  as  follows  at  each  nodal 

points: 
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Here  the  prime  and  the  double  prime  denote  the  first  and  the  second  order  derivative  with 
respect  to  the  spatial  variable  z,  and  the  superscript  ss  indicates  the  steady  state  solution 
which  is  known  as  the  base  flow  for  the  stability  analysis.  The  collection  of  equation  (A.7) 
at  all  nodal  points  along  with  boundary  conditions,  u’*  = w '*+1  = a’  = T’  = 0,  gives  the 
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following  axial  force  balance  in  a matrix  form: 
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By  the  same  way,  the  constitutive  relations  (3.7)  and  (3.8)  are  reduced  in  a matrix 
form  as  follows: 
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By  multiplying  an  inverse  of  H , (i.e.,  H ) on  both  sides  of 

J f J b ==(m— l)x(«— 1)  * v ’ =(n— l)x(«— 1)  ' 

equation  (A. 8),  the  axial  velocity  is  expressed  in  terms  of  other  variables  as  follows: 
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By  substituting  (A.  11)  into  (A. 6),  (A. 9),  and  (A.  10),  the  following  set  of  algebraic 
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equations  are  derived  for  the  eigenvalue  problem: 
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Equations  (A.12)-(A.  14)  can  be  combined  together  to  give 
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Equation  (A.  15)  or  (A.  16)  is  a standard  form  of  an  eigenvalue  problem  for  which  the 
eignevalue  X is  to  be  determined.  In  the  present  study  it  is  solved  by  the  DEN'LRG 

subroutine  in  the  IMSL  package. 


APPENDIX  B 

THE  MAXIMUM  DRAW  RATIO 


The  analytic  expression  for  the  maximum  attainable  draw  ratio  above  which  the 
axial  stress  becomes  infinity  is  discussed  in  this  section.  The  maximun  draw  ratio  for  the 
single-phase  flow  of  the  generalized  UCM  fluid  was  first  obtained  by  Fisher  and  Denn 
(1976)  which  is  described  as 

Dr  <\  + De~Vn  (B.l) 

When  the  power-law  index  n for  the  shear  viscosity  is  set  to  1 , the  result  for  the  single- 
phase UCM  fluid  is  obtained. 

By  a similar  approach,  Park(1990)  obtained  the  analytic  expression  for  the 
maximum  draw  ratio  for  the  two-phase  fiber  spinning  flow: 
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(B.2) 


where  k 3 = 1 ~(t0f)  If  both /and  t0  are  set  to  1,  we  may  find  that  (B.2)  is  reduced  to 

(B.l).  Since  equation  (B.2)  is  implicit  in  Dr,  the  maximum  draw  ratio  is  calculated 
numerically  once  the  values  of  t0,f  (thus  k ) and  De  are  specified.  The  solid  lines  shown 
in  Figure  B are  the  results  of  the  numerical  calculation  for /=0.2,  0.5,  1.0  at  t0  = 10.0  as  a 
function  of  De.  It  should  be  noted  that  the  Dr  value  shown  in  Figure  34  corresponds  to 
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the  infinite  axial  tension.  In  case  of  /=0.2  for  which  the  Dr  vs.  G curve  shows  a 
maximum  at  a finite  value  of  axial  tension,  the  value  of  Dr  which  corresponds  to  the 
maximum  point  is  the  true  maximum  Dr  according  to  the  stability  analysis  of  Chapter  3. 
The  true  maximum  Dr  calculated  from  the  collection  of  numerous  Dr  vs.  G curves  for 
f=  0.2  is  plotted  in  dashed  line  on  Figure  34. 
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Figure  34.  Maximum  attainable  draw  ratio  as  a function  of  Deborah 

number  for/=0.2,  0.5  and  1.0  at  /o=10.0  ( — ),  and  true 
maximum  draw  ratio  as  a function  of  Deborah  number  for 
y=0.2  at  /0=1 0.0  ( ). 


APPENDIX  C 
COMPUTER  PROGRAMS 


The  bicomponent  fiber  spinning  has  been  simulated  by  using  computer  programs. 
These  programs  include  1 main  program(STABCOEX)  and  3 subprgrams  (function 
FUN,  subroutine  RGKG  and  subroutine  DTRIDAG).  There  is  one  additional  subprogram 
which  is  an  IMSL  library  subroutine  (DEVLRG)  called  by  main  program.  The  double 
precision  is  used  throughout  whole  computations.  The  function  of  each  program  is 
described  briefly  as  follows: 

STABCOEX: 

(1)  Solving  inverse  axial  force  G use  shooting  method  as  Dr  is  given. 

(2)  Solving  the  steady-state  solution  use  Runge-Kutta-Gill  method  as  the  parameters 
are  given. 

(3)  Using  linear  stability  analysis  calculate  the  eigenvalues  of  disturbanced  equations 
and  determine  the  stability  status. 

FUN: 

The  functions  of  the  first  order  ODE's  of  equations  (2.37),  (2.38)  and  (2.39). 

RGKG: 

Using  Runge-Kutta-Gill  method  to  solve  the  integration  of  ODE's. 

DTRIDAG: 

Solving  the  inverse  of  a tridigonal  matrix  by  Thomas  algorithm. 

DEVLRG: 

Solving  the  eigenvalues  X of  eigenvalue  problem. 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


C c 
c Program  STABCOEX  c 
c by  c 
c Woei-Shyong  LEE  c 
c c 
c University  of  Florida  c 
c 1993  c 
c c 
c This  program  is  designed  to  solve  the  two-phase  fiber  spinning  c 
c problem,  which  is  assumed  that  the  skin-layer  material  is  an  c 
c upper-convected  Maxwell  fluid  (UCM)  and  the  core-layer  material  c 
c is  a Newtonian  fluid  ( NWF ) . It  is  assumed  that  the  system  is  an  c 
c isothermal,  axisymmetrical , uniaxial  and  neglect  the  inertia,  c 
c surface  tension  and  gravity.  It  can  solve  the  draw  ratio  versus  c 
c inverse  axial  force  distribution,  the  steady-state  solution,  and  c 
c eigenvlues  of  the  disturbanced  equations.  This  program  has  three  c 
c major  parts  which  are:  c 
c (1)  Solving  inverse  axial  force  G use  shooting  method  as  Dr  is  c 
c given.  c 
c (2)  Solving  the  steady-state  solution  use  Runge-Kutta-Gill  c 
c method  as  the  parameters  are  given.  c 
c (3)  Using  linear  stability  analysis  calculate  the  eigenvalues  of  c 
c disturbanced  equations  and  determine  the  stability  status.  c 
c This  computation  is  based  on  the  assumption  of  no  hoop  stress,  c 
c that  is,  E(z)=0  for  all  z c 
c c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


parameter ( ndep=3 , ns=7 , nspace=4200 , n=50 ) 
implicit  double  precision ( a-h, o-z ) 
double  precision  M 

dimension  yy(ns) , y(n+l, 7) ,AMX(n,n-l)  ,BMX(n,n) ,CMX(n, n-1) 
dimension  DMX(n,n) ,EMX(n,n-l) ,FMX(n,n) ,GMX(n,n) , HMX ( n-1 , n-1 ) 
dimension  PMX(n-l,n) ,QMX(n-l,n) , HINV ( n-1 , n-1 ) ,VlMX(n-l,n) 
dimension  V2MX(n-l,n) ,WlMX(n,n) ,W2MX(n,n) ,W3MX(n,n) ,W4MX(n,n) 
dimension  W5MX ( n, n) , W6MX ( n, n) , UMX ( 3*n , 3*n) , tria ( n-1 ) , trib ( n-1 ) 
dimension  trie (n-1) , trif ( n-1 , n-1 ) , RWKSP ( 100000 ) 
complex  eval(3*n) 
common  /worksp/rwksp 
common /block/ f , M, De, G 
data  xO, xmax/0 . 0, 1 . 0/ 
external  DEVLRG,fun 
call  IWKIN( 100000) 
c 

c ndep  - total  number  of  dependent  variables  in  ODEs 

c ns  - total  number  of  variables 

c xO  - initial  value  of  independent  variable 

c h - increment  of  independent  variable 

c xmax  - maximum  of  independent  variable 

c 


open ( 5, file=' STABCOEX. DAT’ , status= • old • ) 
open ( 6, file=' STABCOEX. OUT' , status= ' new ' ) 
nprint  = n 

h = (xmax-xO ) /nspace 


c 

c 

Read 

c 

c 

f 

c 

Dr 

c 

De 

c 

M 

c 

GO 

c 

dGO 

Read  parameters  from  data  file  STABCOEX . DAT  — 


Skin  flow  rate  ratio. 

Draw  ratio. 

Deborah  number. 

Shear  viscosity  ratio  of  core  to  skin. 
Initial  guessed  G value  (G  : inverse  axial 
Initial  setting  of  discritzed  space. 


force ) 


c 

c 

c 
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c 

C = = = = 

c 

c (1) 

c 

c 

c==== 

c 


10 


20 


30 

c 

c==== 

c 

c (2) 
c 
c 
c 

c = = = = 
c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TauzzO  : Initial  extra  shear  stress  at  the  origin  of  the  fiber. 

TO  : Initial  first  normall  stress  at  the  origin  of  the  fiber. 

read( 5, * ) f ,Dr, De , M, GO , dGO , TauzzO , TO 
alpha  = dlog(Dr) 


=================================================================c 

c 

This  part  of  the  program  is  solving  G when  Dr  is  given  by  using  c 
shooting  method.  c 

c 

=============================================== ==================c 


write(*,100)  Dr 
G = GO 
iter  = 0 
iter  = iter  + 1 
x = xO 

yy  ( i ) = i.o 

yy ( 2 ) = yy(l)/3./G/M 
yy(3)  = TauzzO 
do  i = l,nspace 

call  rgkg ( fun, ndep, x , h, yy ) 
x = x0+i*h 
end  do 

ysavenew  = yy(l)  - Dr 
write(*,150)  G,yy(l) 

if ( dabs ( ysavenew) . le . 1 . d-6 ) go  to  30 
if ( iter . eq. 1 ) go  to  20 
ysign  = ysavenew*ysaveold 
if (ysign. It .0.0)  dGO  = dG0*0.1 
if (yy ( 1) .gt .Dr)  then 
G = G + dGO 
else 

if(G.eq.dGO)  dGO  = dG0*0.1 
G = G - dGO 
end  if 

ysaveold  = ysavenew 
go  to  10 

write(*,200)  G,Dr 


========================================================= ========c 

c 

This  part  of  the  program  use  Runge-Kutta-Gill  method  to  find  the  c 
steady-state  solution  of  ODEs  for  two-phase  coextrusion  fiber  c 
spinning.  c 

c 

= = = = = = = = = = =r  = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =:  = =:  = = = = = = = = = = = C 


Solving  steady-state  eqs  and  save  the  solutions  from  yy  vector 
to  y maxtrix  

yy  - vector  of  dependent  variables 

yy(l)  = W(ss),  yy (2 ) = W'(ss),  yy(3)  = Tau/zz(ss) 

yy(4)  = W"(ss),  yy(5)  = Tau' ,zz(ss) 

yy ( 6 ) = t ( ss ) , yy(7)  = T'(ss) 

: node  number 
= 1 : W( ss ) 

2 : W’ (ss) 

3 : Tau,zz(ss) 

4 : W" (ss) 

5 : Tau ,zz' ( ss ) 

6 : T ( ss ) 


y(i,k) 


1 

k 


; ss: steady  state 


**  e» 


c 

c 


7 : T‘ (ss) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c: 

c 


x = xO 

yy ( i ) = l.o 

yy ( 2 ) = yy(l)/3./G/M 
yy(3)  = TauzzO 

yy  ( 4 ) = l./(3.*(l.-f)*G*M*De)  - yy(2 ) **2/yy ( 1)  + yy(2) 

* (2./(3.*(l.-f)*G*M)  - f/(l.-f)/M/De/yy(l) 

- f*yy(3)/(l.-f )/G/M/yy(l)  - l./De/yy(l)) 
yy ( 5 ) = ( 2 . *yy ( 2 ) * (G+De*yy ( 3 ) ) -yy ( 3 ) ) /De/yy ( 1 ) 

yy ( 6 ) = yy(l) 

yy(7)  = (yy(2)*(3.*(G+De*yy(3) )-De*yy(6) )-yy(6) )/De/yy(l) 
k = 1 

do  j = l,ns 

y(*,j)  = yy(j) 

end  do 

> i = l,nspace 

iprint=nspace/nprint *int ( i*nprint/nspace) 
call  rgkg( fun, ndep, x, h, yy ) 

yy ( 4 ) = 1 • / (3.*(l.-f)*G*M*De)  - yy ( 2 ) * *2/yy ( 1 ) + yy(2) 

* (2./(3.*(l.-f )*G*M)  - f/(l.-f )/M/De/yy(l) 

- f*yy(3)/(l.-f)/G/M/yy(l)  - l./De/yy(l)) 

YY (5 ) = ( 2 . *yy ( 2 ) * (G+De*yy ( 3 ) ) -yy ( 3 ) ) /De/ yy ( 1 ) 
yy ( 6 ) = yy(l)/f  - 3 . * ( 1 . -f ) *G*M*yy (2 ) /f 

yy(7)  = (yy(2)*(3.*(G+De*yy(3) )-De*yy(6) )-yy(6) )/De/yy(l) 
if ( i . eq. iprint ) then 
k = k + 1 
do  j = l,ns 

y(k,j)  = yy ( j ) 
end  do 
end  if 
x = x0+i*h 
end  do 

write ( * , 250 ) f , G, TauzzO , M, Dr , TO, De, alpha, nspace , n 
write ( 6 , 250 ) f , G, TauzzO, M, Dr , TO, De, alpha, nspace, n 
do  k = l,n+l 

x = xO  + (k-1) /float (n) 
write ( * , 300)  x,y(k,l),y(k,3),y(k,6) 
write (6, 300)  x,y(k,l),y(k,3),y(k,6) 
end  do 
write ( * , 350 ) 
write ( 6 , 350 ) 
do  k = l,n+l 

x = xO  + ( k-1 ) /float ( n ) 

wr ite ( * , 400 ) x,y(k,2),y(k,4),y(k,5),y(k,7) 
write (6, 400)  x,y(k,2),y(k,4),y(k,5),y(k,7) 
end  do 
write( * , 450 ) 


c 

(3)  This  part  of  the  program  is  used  to  solve  all  the  eigenvalues  of  c 
transition  equations  of  two  phase  coextrusion  fiber  spinning  c 

(skin:  UCM,  core:NWF).  The  first  four  largest  eigenvalues  of  c 

disturbanced  equations  are  printed  out  to  determine  whether  the  c 
system  is  stable  or  not.  The  system  is  stable  if  all  real  parts  c 
of  eigenvalues  are  less  than  zero  while  the  system  is  unstable  c 
if  at  least  one  real  part  of  eigenvalue  is  greater  than  zero.  If  c 
the  real  part  of  eigenvalue  is  zero,  it  is  called  neutral  stable. c 
An  IMSL  subroutine  DVLRG  is  called  in  this  program.  c 

c 


c 


y(i,k)  : i : node  number 


; ss: steady  state 


k = 1 
2 

3 

4 

5 

6 
7 

h=l . O/float ( n ) 


W(  ss ) 

W ' ( ss ) 

Tau , zz ( ss ) 
W" (88) 

Tau , zz ' ( ss ) 
T ( ss ) 

T' (ss) 


AMX (n, n-1) *w(n-l, l)+BMX(n,n) *a(n, 1 

AMX(n, n-1) 

do  i = 1 , n 
do  j = l,n-l 
AMX(i,j)  =0.0 
end  do 
end  do 
do  i = 2,n+l 

if ( i.ge. 3 . and. i. le. (n-1 ) ) then 
AMX(i-l,i-2)  = 0 . 5*y ( i , 1 ) 

AMX(i-l,i-l)  = 0.0 

AMX ( i-1 , i ) = -0 . 5*y ( i , 1 ) 

else 

if(i.eq.2)  then 

AMX ( i-1 , i-1 ) = 0.0 

AMX(i-l,i)  = -0.5*y(i,l) 
end  if 

if(i.eq.n)  then 

AMX ( i-1 / i-2 ) = 0 . 5*y ( i , 1 ) 

AMX ( i-1 , i-1 ) = 0.0 

end  if 

if ( i . eq. ( n+1 ) ) then 
AMX ( i-1 , i-2 ) = Dr 
end  if 
end  if 
end  do 


BMX ( n,  n ) 

do  i = l,n 
do  j = l,n 

BMX ( i , j ) = 0.0 
end  do 
end  do 
do  i = 2, n+1 

if ( i . ge . 3 . and . i . le . n ) then 

BMX ( i-1 , i-2 ) = 0 . 5*y ( i , 1 ) 

BMX ( i-1 , i-1 ) = 0.0 

BMX ( i-1 , i ) = -0 . 5*y ( i , 1 ) 

else 

if(i.eq.2)  then 

BMX ( i-1 , i-1 ) = 0.0 

BMX ( i-1 , i ) = -0 . 5*y ( i , 1 ) 

end  if 

if ( i . eq. ( n+1 ) ) then 
BMX ( i-1 f i-2 ) = Dr 
BMX ( i-1 , i-1 ) = -Dr 
end  if 
end  if 
end  do 


= lambda* h* a ( n, 1 ) 


CMX(n, n-1) *w(n-l, l)+DMX(n/n) *Tau(n, 1 ) =lambda*h*Tau ( n,  1) 


CMX ( n, n-1 ) 


do  i = 1 , n 
do  j = l,n-l 
CMX ( i , j ) =0.0 
end  do 
end  do 

do  i = 2, n+1 

aelm  = -y ( i, 1 ) * ( 1 . +G/De/y ( i, 3 ) ) 
belm  = 2 . *h*y ( i, 2 ) +2 . *h*G*y ( i, 2 ) /De/y ( i,  3 ) 

& -h*y(i,l)*y(i,5)/y(i,3) 

if ( i.ge. 3 . and. i. le. (n-1 ) ) then 


CMX ( i-1 , i-2 ) = aelm 
CMX ( i-1 , i-1 ) = belm 
CMX ( i-1 , i ) = -aelm 

else 

if(i.eq.2)  then 

CMX ( i-1 , i-1 ) = belm 
CMX ( i-1 , i ) = -aelm 

end  if 

if(i.eq.n)  then  * 

CMX ( i-1 , i-2 ) = aelm 
CMX (i-1, i-1)  = belm 
end  if 


if ( i . eq. ( n+1 ) ) CMX(i-l,i-2)  = 2.*aelm 
end  if 
end  do 


■ DMX ( n, n) 

do  i = l,n 
do  j = l,n 

DMX ( i , j ) = 0.0 
end  do 
end  do 

do  i = 2, n+1 

celm  = 0 . 5*y ( i, 1 ) 

delm  = 2 . *h*y ( i, 2 ) -h/De-h*y ( i,l)*y(i/5)/y(i/3) 
if ( i. ge. 3 . and. i. le. n)  then 
DMX ( i-1 , i-2 ) = celm 
DMX ( i-1 , i-1 ) = delm 
DMX ( i-1 , i ) = -celm 

else 

if(i.eq.2)  then 

DMX ( i-1 , i-1 ) = delm 
DMX ( i-1 , i ) = -celm 

end  if 

if ( i. eq. ( n+1 ) ) then 

DMX ( i-1 , i-2 ) = 2 . *celm 
DMX(i-l,i-l)  = delm-2.*celm 
end  if 
end  if 
end  do 

EMX ( n , n-1 ) *w ( n-1 , 1 ) +h*FMX ( n , n ) *Tau ( n, 1 ) 

+GMX (n,n)*T(n,l)  = lambda*h*T ( n, 1 ) 

EMX ( n, n-1 ) 


do  i = l,n 
do  j = l,n-l 
EMX ( i , j ) =0.0 
end  do 


e> 


end  do 
do  i = l,n 

eelm  = 0 . 5*y ( i+1 , 1 ) * ( 1 . -3 . *y ( i+1 , 3 ) /y ( i+1 , 6 ) 

& -3.*G/De/y(i+l,6) ) 

felm  = 3 . *h*y ( i+1 , 2 ) *y ( i+1 , 3 ) /y ( i+1 , 6 ) 

+3. *h*G*y(i+l,2) /De/y(i+l, 6) 

-h*y (i+1, 1) *y (i+1, 7 ) /y ( i+1 , 6 ) -h*y ( i+1 , 2 ) 
if ( i . ge . 2 . and . i . le . ( n-2 ) ) then 
EMX (i, i-1)  = eelm 
EMX ( i , i ) = felm 

EMX ( i , i+1 ) = -eelm 
else 

if(i.eq.l)  then 

EMX(i,i)  = felm 
EMX ( i , i+1 ) = -eelm 
end  if 

if ( i . eq. ( n-1 ) ) then 
EMX(i, i-1)  = eelm 
EMX ( i , i ) = felm 

end  if 

if(i.eq.n)  EMX(i,i-l)  = 2.*eelm 
end  if 
end  do 


FMX ( n , n ) 

do  i = l,n 
do  j = l,n 

FMX ( i , j ) =0.0 

if(i.eq.j)  FMX ( i , j ) = 3. *y (i+l# 3) *y ( i+lf 2 ) /y ( i+l# 6) 
end  do 
end  do 

GMX ( n , n ) 

do  i = l,n 
do  j = 1 , n 

GMX ( i , j ) =0.0 
end  do 
end  do 
do  i = l,n 

gelm  = 0 . 5*y ( i+1, 1 ) 

pelm  = -h/De-h*y ( i+1 ,l)*y(i+l,7)/y(i+l,6)-h*y(i+l,2) 
if ( i . ge. 2 . and. i . le. ( n-1 ) ) then 
GMX(i/i-l)  = gelm 
GMX(i,i)  = pelm 
GMX (i, i+1)  = -gelm 
else 

if(i.eq.l)  then 

GMX(i,i)  = pelm 
GMX (i, i+1)  = -gelm 
end  if 

if(i.eq.n)  then 

GMX ( i , i-1 ) = 2*gelm 
GMX(i/i)  = pelm-2*gelm 
end  if 
end  if 
end  do 

HMX ( n-1 / n-1 ) *w( n-1 f 1 ) =h*PMX (n-l,n)*a(n/l) +h*QMX (n-lfn)*T(n/l) 
HMX ( n-1 , n-1 ) 


do  j = 1 , n- 1 
HMX(i,j)  =0.0 
end  do 
end  do 

do  i = 1 , n-1 

qelml  = 3 . * ( 1 . -f ) *G*M* ( y ( i+1 , 1 ) * *2-0 . 5*h*y ( i+1 , 1 ) *y ( i+1 , 2 ) ) 
relm  = 3 . * ( 1 . -f ) *G*M* ( -2 . *y ( i+1 , 1 ) * *2+h* *2*y ( i+1 , 1 ) *y ( i+1 , 4 ) 
Sc  -h**2*y ( i+1  / 2 ) **2  ) 

qelm2  = 3. * ( l.-f ) *G*M* (y ( i+1, 1 ) **2+0. 5*h*y ( i+1 , 1 ) *y ( i+1, 2 ) ) 
if ( i . ge. 2 . and. i . le. ( n-2 ) ) then 
HMX ( i , i-1 ) = qelml 
HMX ( i , i ) = relm 

HMX (i, i+1)  = qelm2 
else 

if(i.eq.l)  then 

HMX ( i , i ) = relm 

HMX (i, i+1)  = qelm2 
end  if 

if ( i . eq. ( n-1 ) ) then 
HMX ( i , i-1 ) = qelml 
HMX(i,i)  = relm 
end  if 
end  if 
end  do 

PMX(n-l,n) 

do  i = l,n-l 
do  j = 1 , n 

PMX(i,j)  =0.0 
end  do 
end  do 

do  i = l,n-l 

selm  = 0.5*y(i+l,l)*(f*y(i+l,6)+3.*(l.-f)*G*M*y(i+l,2) ) 
telm  = h* ( f * (y ( i+1, 6) *y ( i+1, 2 ) -y ( i+1, 1) *y ( i+1, 7 ) ) 

Sc  +3.*(l.-f)*G*M*(y(i+l,2)**2-y(i+l,l)*y(i+l,4)  ) ) 

if(i.eq.l)  then 


PMX ( i f i ) 

= telm 

PMX ( i , i+1 ) 

= -selm 

else 

PMX (i, i-1) 

= selm 

PMX ( i , i ) 

= telm 

PMX (i, i+1) 

= -selm 

end  if 

end  do 

QMX ( n-1, n) 

do  i = l,n-l 

do  j = 1 r n 

QMX(i,j)  = 

o 

• 

o 

end  do 


end  do 

do  i = l,n-l 

uelm  = 0. 5*f *y (i+1, 1 ) *y ( i+1, 6) 

velm  = h*f * (y ( i+1 , 6 ) *y ( i+1 , 2 ) -y ( i+1 , 1 ) *y ( i+1 , 7 ) ) 
if(i.eq.l)  then 

QMX ( i , i ) = velm 

QMX(i,i+l)  = -uelm 
else 

QMX (1,1-1)  = uelm 
QMX ( i , i ) = velm 

QMX (i, i+1)  = -uelm 
end  if 


end  do 


w ( n-1 , 1 ) = h*HINV(n-l,n-l)  *PMX  ( n-1 , n ) *a ( n , 1) 

+ h*HINV(n-l,n-l) *QMX(n-l,n) *T(n, 1 ) 

Compute  HINV ( n-1 , n-1 ) , inverse  of  HMX  matrix  by  calling 
subroutine  DTRIDAG 

do  i = l,n-l 

if ( i . ge. 2 . and. i . le. (n-2  ) ) then 
tria(i)  = HMX(i,i-l) 
trib(i)  = HMX ( i , i ) 
trie ( i ) = HMX ( i , i+1 ) 
else 

if(i.eq.l)  then 

trib(i)  = HMX ( i / i ) 
trie ( i ) = HMX ( i , i+1 ) 
end  if 

if ( i . eq. (n-1 ) ) then 
tria(i)  = HMX(i,i-l) 
trib(i)  = HMX(i/i) 
end  if 
end  if 

do  j = l,n-l 

trif ( i , j ) = 0.0 
if(i.eq.j)  trif(i,j)  = 1.0 
end  do 
end  do 

call  DTRIDAG(tria, trib, trie, trif , n-1 , HINV ) 

HINV ( n-1 , n-1 ) *PMX ( n-1 , n)  — > VlMX(n-l,n) 

do  i = l,n-l 
do  j = l,n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + HINV( i, k) *PMX ( k, j ) 
end  do 

VlMX(irj)  = sum 
end  do 
end  do 

HINV ( n-1 , n-1 ) *QMX(n-l, n)  — > V2MX(n-l/n) 

do  i = lfn-l 
do  j = l,n 
sum  = 0.0 
do  k = 1,0-1 

sum  = sum  + HINV( i, k) *QMX ( k, j ) 
end  do 

V2MX(i,j)  = sum 
end  do 
end  do 

h*AMX ( n, n-1 ) *V1MX ( n-1 , n ) + BMX(n,n)  — > WlMX(n,n) 

do  i = l,n 
do  j = 1 , n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + AMX  ( i , k)  *V1MX ( k,  j ) 
end  do 

W1MX ( i , j ) = h*sum  + BMX(i,j) 
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end  do 
end  do 
c 

c h*AMX(n,n-l)*V2MX(n-l,n)  — > W2MX(n,n) 

c 

do  i = l,n 
do  j = l,n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  4-  AMX  ( i , k)  *V2MX  ( k,  j ) 
end  do 

W2MX(i,j)  = h*sum 
end  do 
end  do 
c 

c h*CMX ( n, n-1 ) *V1MX ( n-1 , n ) — > W3MX(n,n) 

c 

do  i = l,n 
do  j = l,n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + CMX ( i , k) *V1MX ( k, j ) 
end  do 

W3MX(i,j)  = h*sum 
end  do 
end  do 
c 

c h*CMX(n,n-l)*V2MX(n-l,n)  — > W4MX(n,n) 

c 

do  i = l,n 
do  j = l,n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + CMX ( i , k) *V2MX ( k, j ) 
end  do 

W4MX(i,j)  = h*sum 
end  do 
end  do 
c 

c h*EMX(n,n-l)*VlMX(n-l,n)  — > W5MX(n/n) 

c 

do  i = l,n 
do  j = l,n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + EMX ( i , k) *V1MX ( k, j ) 
end  do 

W5MX(i,j)  = h*sum 
end  do 
end  do 
c 

c h*EMX ( n , n-1 ) *V2MX ( n-1 , n ) + GMX(n,n)  — > W6MX(n,n) 

c 

do  i = l,n 
do  j = 1 , n 
sum  = 0.0 
do  k = 1 , n-1 

sum  = sum  + EMX ( i , k) *V2MX ( k, j ) 
end  do 

W6MX(i,j)  = h*sum  + GMX(i,j) 
end  do 
end  do 


c 
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C 

C 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

100 

150 

200 

250 


300 

350 

400 

450 


UMX ( 3n, 3n ) = 


W1MX ( n , n ) 
{ W3MX(n,n) 
W5MX ( n , n ) 


OMX ( n, n ) 
DMX ( n , n ) 
h*FMX(n, n) 


do  i = l,n 
do  j = l,n 

W1MX( i, j ) 

0.0 

W2MX(i, j) 
W3MX(i, j ) 
DMX ( i , j ) 
W4MX(i, j) 
W5MX(i, j) 
h*FMX(i, j ) 
W6MX ( i , j ) 

end  do 
end  do 


UMX ( i,  j ) 

UMX ( i, j+n ) 

UMX ( i, j+2  *n ) 

UMX ( i+n, j ) 

UMX ( i+n , j +n ) 

UMX ( i+n, j+2*n) 

UMX( i+2*n, j ) 
UMX(i+2*n, j+n) 

UMX (i+2*n,j+2*n)  = 


W2MX ( n , n) 
W4MX ( n , n ) } 

W6MX(n/n) 


- Call  subroutine  DEVLRG  to  compute  eigenvalues  lambda 
of  UMX(3n,3n)  * x(3n,l)  = lambda  * h * x(3n,l) 

call  DEVLRG ( 3*n, UMX, 3*n, eval ) 

write ( * , 500 ) 
write (6, 500 ) 
do  j = 3*n,3*n-6,-2 

evr  = dble(eval (2* j-1 ) ) /h 
evi  = dble ( eval ( 2* j ) ) /h 
write(*,550)  ( 3*n- j+2 ) /2 , evr , evi 

write(6,550)  ( 3*n- j +2 ) /2 , evr , evi 

if(j.eq.3*n)  evalue  = evr 
end  do 

if (evalue. It . 0 . 0 ) then 
write ( * , 600 ) 
write (6, 600 ) 
end  if 

if ( evalue. gt . 0 . 0 ) then 
write ( * , 650 ) 
write (6,  650 ) 
end  if 

if (dabs (evalue) . le. l.d-6)  then 
write ( *,  700 ) 
write (6,  700 ) 
end  if 
write ( * , 750 ) 

format ( /lx, ' Program  is  looking  for  a G value  when  Dr  =',fl2.6 
& //17x , * G ' , 18x,  'Dr  1 /) 

format ( lx, 2f 20 . 6 ) 

format (/'  Congrulat ions ! , You  finally  get  the  solution  which  is 
& / ' G =',fl2. 6,'  when  Dr  =',fl2.6 

& /’  Now,  program  is  solving  the  steady-state  solution.'/) 

format ( //3x,  ' f = ' , f 12 • 6, 5x, ' G =',fl2.6,5x,’ Tau, zz (0)  = ' , f 12. 6 

& /3x , ' M *' f fl2.6,5x, 'Dr  = ' , f 12 . 6, 5x, ' T ( 0 ) =',fl2.6 

Sc  /3x,'De  = ' , f 12 . 6 , 5x,  ' alpha  = ' , f 12 . 6, 5x,  ' nsapce  =',i5 

Sc  /3x,'n  = ' , i5 

Si  ///3x,'The  steady-state  solution  is:' 

Si  //10x/'z',llx,'w'/20x,'  Tau,  zz'^lOx,  'T'/) 

format (2x, 2 f 12 .6,12x,2fl2.6) 

format ( //10x,  ' z ' , lOx, ' w’  ' ' , lOx, ' w" ' , 7x,  ' Tau, zz ' ' ' , 9x,  ' T ' ' ' / ) 
format ( 2x,  5f 12 . 6 ) 

format(/lx,'  Now,  program  is  computing  the  eigenvalues.' 


• • 
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& /lx,  * Wait  a moment  Be  patient  ! ' ) 

500  format ( //lx , ' The  eigenvalues  are:' 

& //29x , ' 1 , real ' , 14x , ' 1 , imag ' / ) 

550  format ( lOx, i4,  ' :',2f20.6) 

600  format ( /3x ,' Stable : all  eigenvalues  are  less  than  zero.'/) 

650  format ( /3x ,' Unstable:  at  least  one  eigenvalue  is  greater  ' 

& 'than  zero'/) 

700  format ( /3x ,' Neutral  Stable:  the  eigenvalue  is  close  to  zero.'/) 

750  format(/lx,'  Computation  is  completed.  If  you  want  to  review  the' 

& /lx,'  solution,  please  input  TYPE  STABCOEX . OUT ' / ) 

end 


c 

c 

c 

c 

c 

c 

c 

c 

c 


1 

2 


c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c — 


c 

c 


function  fun(j,n,x,y) 

j - sequential  number  of  dependent  variables 

n - total  number  of  dependent  variables 

x - independent  variable 

y - vector  dependent  variables 

implicit  double  precision ( a-h, o-z ) 

double  precision  M 

dimension  y(n) 

common/block/f ,M, De,G 

go  to  (1,2,3) , j 

fun  = y ( 2 ) 

return 

fun  = l./(3.*(l.-f)*G*M*De)  “ y(2)**2/y(l)  + y(2) 
& * (2./(3.*(l.-f)*G*M)  - f/(l.-f )/M/De/y(l) 

& - f *y ( 3 ) / (l.-f)/G/M/y(l)  - l./De/y(l)) 

return 

fun  = ( 2 . *y ( 2 ) * (G+De*y ( 3 ) ) -y ( 3 ) ) /De/y ( 1 ) 

return 

end 


c 

c 


Subroutine  RGKG ( f , n, x, h, y ) 

This  subprogram  use  Runge-Kutta-Gill  method  to  solve  the 
system  of  simultaneous  first  order  ordinary  differential 
equations . 


n 

x 

h 

y 


- total  number  of  dependent  variables 

- independent  variable 

- increment  of  independent  variable 

- vector  of  dependent  variables 


c — 


implicit  double  precision ( a-h, o-z ) 

dimension  y ( n ) , ybar ( 50 ) , phi ( 50 ) 

double  precision  k(50,4) 

a=dsqrt ( 2 . dO ) 

a=l .0/a 

b=l . 0-a 

c=l . 0+a 

step  1 

do  j=l,n 

k( j,l)=f ( j,n,x,y) 
ybar ( j ) =y ( j ) 
phi( j ) =k( j , 1 ) 
end  do 

step  2 

do  j=l,n 
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y ( j )=ybar ( j )+0. 5*h*k( j ,1) 
end  do 
x=x+0 . 5 *h 
do  j = l,n 

k( j,2)=f ( j,n,x,y) 

phi ( j ) =phi (j)+2.0*b*k(j,2) 
end  do 

c — step  3 

do  j = l,n 

y ( j ) =ybar ( j ) +b*h* ( a*k ( j , 1 ) +k ( j , 2 ) ) 
end  do 
do  j=l,n 

k( j,3)=f ( j,n,x,y) 

phi ( j ) =phi ( j )+2.0*c*k( j,3) 
end  do 

c — step  4 

do  j=l,n 

y ( j ) =ybar ( j ) -a*h*k( j , 2 ) +c*h*k ( j , 3 ) 
end  do 
x=x+0 . 5*h 
do  j = l,n 

k(j/4)=f(j,n/x/y) 
end  do 

c — step  5 

do  j = l,n 

y ( j ) =ybar ( j ) +h* ( phi (j)+k(j,4))/6.0 
end  do 
return 
end 

c 


c c 

c c 

Subroutine  DTRIDAG(a,b, c, f ,n,x) 
c 

c This  subprogram  use  Thomas  Algorithm  to  calculate  the  inverse 

c matrix  of  a tridigonal  matrix, 

c 

c n total  magnitude  of  x matrix 

c x - inverse  matrix 

c f identity  matrix 

c a,b,c  - tridigonal  elements 

c 

implicit  double  precision ( a-h, o-z ) 


dimension  a(n),b(n),c(n),f(n,n),x(n,n) 
do  i = 2 , n 

b ( i ) = b(i)  - c(i-l)*a(i)/b(i-l) 
do  j = 1 , n 

f(i,j)  = f(i,j)  - f(i-l,j)*a(i)/b(i-l) 
end  do 
end  do 
do  j = l,n 

x(n,j)  = f(n,j)/b(n) 
do  i = l,n-l 

x(n-i, j ) = (f (n-i, j )-c(n-i)*x(n-i+l, j ) ) /b(n-i) 
end  do 
end  do 
return 
end 


Example  1: 
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f 

M 

De 

n 


0.200000 

G 

0.106984 

Tau,zz(0)  = 10 

1.000000 

Dr 

20.000000 

T ( 0 ) = 1 

0.020000 

alpha  = 

2.995732 

nsapce  = 4200 

50 

The  steady-state  solution  is: 


z w 


0.000000 

1.000000 

0.020000 

1.057770 

0.040000 

1.117991 

0.060000 

1.184732 

0.080000 

1.257842 

0.100000 

1.336898 

0.120000 

1.421722 

0.140000 

1.512359 

0.160000 

1.608999 

0.180000 

1.711926 

0.200000 

1.821479 

0.220000 

1.938046 

0.240000 

2.062047 

0.260000 

2.193935 

0.280000 

2.334193 

0.300000 

2.483334 

0.320000 

2.641904 

0.340000 

2.810478 

0.360000 

2.989664 

0.380000 

3.180108 

0.400000 

3.382485 

0.420000 

3.597513 

0.440000 

3.825943 

0.460000 

4.068568 

0.480000 

4.326220 

0.500000 

4.599773 

0.520000 

4.890142 

0.540000 

5.198287 

0.560000 

5.525210 

0.580000 

5.871959 

0.600000 

6.239625 

0.620000 

6.629343 

0.640000 

7.042294 

0.660000 

7.479700 

0.680000 

7.942827 

0.700000 

8.432981 

0.720000 

8.951507 

0.740000 

9.499783 

0.760000 

10.079224 

0.780000 

10.691272 

0.800000 

11.337393 

0.820000 

12.019073 

0.840000 

12.737812 

0.860000 

13.495114 

0.880000 

14.292486 

0.900000 

15.131421 

0.920000 

16.013397 

0.940000 

16.939861 

0.960000 

17.912222 

0.980000 

18.931837 

1.000000 

19.999999 

Tau , zz 

T 

10.000000 

1.000000 

4.634612 

1.601969 

2.472781 

1.519806 

1.603187 

1.429720 

1.259816 

1.402097 

1.139957 

1.424216 

1.120466 

1.478581 

1.149465 

1.553714 

1.204371 

1.643254 

1.275113 

1.743965 

1.357177 

1.854349 

1.448618 

1.973856 

1.548738 

2.102465 

1.657478 

2.240459 

1.775130 

2.388314 

1.902202 

2.546640 

2.039351 

2.716153 

2.187349 

2.897661 

2.347073 

3.092063 

2.519502 

3.300346 

2.705714 

3.523593 

2.906896 

3.762987 

3.124348 

4.019822 

3.359497 

4.295509 

3.613901 

4.591590 

3.889265 

4.909748 

4.187454 

5.251819 

4.510505 

5.619808 

4.860643 

6.015903 

5.240298 

6.442491 

5.652120 

6.902177 

6.098998 

7.397801 

6.584083 

7.932456 

7.110801 

8.509510 

7.682881 

9.132631 

8.304374 

9.805800 

8.979674 

10.533344 

9.713544 

11.319951 

10.511135 

12.170697 

11.378013 

13.091071 

12.320178 

14.086993 

13.344086 

15.164837 

14.456670 

16.331454 

15.665354 

17.594185 

16.978072 

18.960880 

18.403275 

20.439907 

19.949942 

22.040164 

21.627579 

23.771073 

23.446215 

25.642585 

25.416391 

27.665161 

27.549140 

29.849761 

000000 

000000 
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z 

0.000000 

0.020000 

0.040000 

0.060000 

0.080000 

0.100000 

0.120000 

0.140000 

0.160000 

0.180000 

0.200000 

0.220000 

0.240000 

0.260000 

0.280000 

0.300000 

0.320000 

0.340000 

0.360000 

0.380000 

0.400000 

0.420000 

0.440000 

0.460000 

0.480000 

0.500000 

0.520000 

0.540000 

0.560000 

0.580000 

0.600000 

0.620000 

0.640000 

0.660000 

0.680000 

0.700000 

0.720000 

0.740000 

0.760000 

0.780000 

0.800000 

0.820000 

0.840000 

0.860000 

0.880000 

0.900000 

0.920000 

0.940000 

0.960000 

0.980000 

1.000000 


w' 

3.115737 

2.871840 

3.170380 

3.500486 

3.806741 

4.097408 

4.385423 

4.679900 

4.986537 

5.308954 

5.649649 

6.010550 

6.393315 

6.799487 

7.230575 

7.688107 

8.173644 

8.688800 

9.235248 

9.814723 

10.429024 

11.080013 

11.769617 

12.499819 

13.272661 

14.090235 

14.954677 

15.868161 

16.832888 

17.851076 

18.924947 

20.056714 

21.248561 

22.502628 

23.820987 

25.205621 

26.658398 

28.181040 

29.775097 

31.441911 

33.182583 

34.997934 

36.888467 

38.854327 

40.895261 

43.010572 

45.199084 

47.459094 

49.788336 

52.183948 

54.642429 


w" 

-58.246921 

10.215344 

16.815563 

15.921666 

14.804915 

14.373479 

14.502957 

14.991043 

15.702069 

16.560006 

17.525314 

18.578514 

19.710779 

20.918903 

22.202669 

23.563482 

25.003633 

26.525907 

28.133352 

29.829141 

31.616480 

33.498534 

35.478368 

37.558885 

39.742770 

42.032424 

44.429902 

46.936832 

49.554344 

52.282979 

55.122600 

58.072291 

61.130257 

64.293707 

67.558748 

70.920260 

74.371781 

77.905382 

81.511554 

85.179090 

88.894980 

92.644316 

96.410210 

100.173730 

103.913866 

107.607513 

111.229496 

114.752627 

118.147810 

121.384187 

124.429337 


1 , real 

-0.293327 

-3.302049 

-4.336087 

-4.994560 


Tau , zz ' 

-404.351917 

-164.862637 

-66.227565 

-26.576420 

-10.075412 

-2.857730 

0.507147 

2.216964 

3.194899 

3.843956 

4.347221 

4.791534 

5.220214 

5.656306 

6.113180 

6.599523 

7.121773 

7.685345 

8.295277 

8.956590 

9.674492 

10.454514 

11.302608 

12.225219 

13.229345 

14.322602 

15.513275 

16.810375 

18.223693 

19.763855 

21.442373 

23.271702 

25.265282 

27.437591 

29.804179 

32.381707 

35.187967 

38.241895 

41.563573 

45.174204 

49.096077 

53.352504 

57.967726 

62.966795 

68.375416 

74.219755 

80.526198 

87.321078 

94.630342 

102.479183 

110.891612 


T * 

90.356387 

1.244689 

-5.736008 

-2.937887 

0.027080 

2.034294 

3.308141 

4.153921 

4.774288 

5.284948 

5.749155 

6.201559 

6.661776 

7.141637 

7.648975 

8.189614 

8.768423 

9.389900 

10.058495 

10.778807 

11.555717 

12.394474 

13.300771 

14.280804 

15.341333 

16.489731 

17.734048 

19.083057 

20.546315 

22.134219 

23.858056 

25.730060 

27.763463 

29.972542 

32.372661 

34.980306 

37.813108 

40.889864 

44.230527 

47.856198 

51.789078 

56.052415 

60.670405 

65.668079 

71.071141 

76.905776 

83.198414 

89.975448 

97.262908 

105.086087 

113.469116 


The  eigenvalues  are: 


1 

2 

3 

4 


1 , imag 

13.731118 

33.178902 

61.466445 

69.644864 


Stable:  all  eigenvalues  are  less  than  zero. 


Example  2: 
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f 

M 

De 

n 


0.200000 

G 

0.091367 

Tau,zz(0)  = 10.000000 

1.000000 

Dr 

30.000000 

T ( 0 ) = 1.000000 

0.020000 

0 

alpha  = 

3.401197 

nsapce  = 4200 

The  steady-state  solution  is: 


z 

w 

Tau , zz 

T 

0.000000 

1.000000 

10.000000 

1.000000 

0.020000 

1.066778 

4.721847 

1.699080 

0.040000 

1.136793 

2.568179 

1.611991 

0.060000 

1.215584 

1.694314 

1.514432 

0.080000 

1.303112 

1.348332 

1.486551 

0.100000 

1.398962 

1.231238 

1.515818 

0.120000 

1.503018 

1.220327 

1.583662 

0.140000 

1.615454 

1.263287 

1.677692 

0.160000 

1.736643 

1.337071 

1.791006 

0.180000 

1.867088 

1.431344 

1.920151 

0.200000 

2.007387 

1.541555 

2.063661 

0.220000 

2.158212 

1.665921 

2.221204 

0.240000 

2.320296 

1.804068 

2.393120 

0.260000 

2.494435 

1.956397 

2.580175 

0.280000 

2.681480 

2.123782 

2.783435 

0.300000 

2.882342 

2.307430 

3.004206 

0.320000 

3.097993 

2.508810 

3.244007 

0.340000 

3.329466 

2.729628 

3.504565 

0.360000 

3.577857 

2.971825 

3.787817 

0.380000 

3.844328 

3.237580 

4.095928 

0.400000 

4.130109 

3.529329 

4.431305 

0.420000 

4.436498 

3.849787 

4.796625 

0.440000 

4.764865 

4.201968 

5.194854 

0.460000 

5.116649 

4.589217 

5.629284 

0.480000 

5.493361 

5.015244 

6.103558 

0.500000 

5.896583 

5.484153 

6.621710 

0.520000 

6.327968 

6.000479 

7.188198 

0.540000 

6.789236 

6.569233 

7.807948 

0.560000 

7.282174 

7.195935 

8.486388 

0.580000 

7.808630 

7.886658 

9.229499 

0.600000 

8.370509 

8.648076 

10.043849 

0.620000 

8.969766 

9.487498 

10.936645 

0.640000 

9.608399 

10.412917 

11.915772 

0.660000 

10.288440 

11.433049 

12.989833 

0.680000 

11.011943 

12.557368 

14.168190 

0.700000 

11.780972 

13.796140 

15.460997 

0.720000 

12.597589 

15.160453 

16.879226 

0.740000 

13.463830 

16.662231 

18.434686 

0.760000 

14.381698 

18.314245 

20.140038 

0.780000 

15.353133 

20.130107 

22.008784 

0.800000 

16.379996 

22.124256 

24.055255 

0.820000 

17.464047 

24.311916 

26.294571 

0.840000 

18.606914 

26.709044 

28.742590 

0.860000 

19.810074 

29.332251 

31.415828 

0.880000 

21.074823 

32.198700 

34.331359 

0.900000 

22.402251 

35.325978 

37.506687 

0.920000 

23.793215 

38.731944 

40.959592 

0.940000 

25.248313 

42.434544 

44.707950 

0.960000 

26.767859 

46.451604 

48.769522 

0.980000 

28.351861 

50.800595 

53.161722 

1.000000 

30.000001 

55.498379 

57.901360 
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z 

w ' 

w" 

Tau , zz ' 

0.000000 

3.648298 

-79.860456 

-393.700713 

105 

0.020000 

3.315218 

13.570954 

-163.571379 

1 

0.040000 

3.713942 

22.611396 

-66.326709 

-6 

0.060000 

4.162238 

21.851353 

-26.803910 

-3 

0.080000 

4.586831 

20.728525 

-10.082836 

0 

0.100000 

4.997247 

20.455299 

-2.571928 

2, 

0.120000 

5.409904 

20.910145 

1.075120 

4, 

0.140000 

5.836893 

21.853995 

3.041136 

5, 

0.160000 

6.286208 

23.121923 

4.256333 

6, 

0.180000 

6.763296 

24.620213 

5.135313 

6, 

0.200000 

7.272221 

26.299990 

5.871948 

7, 

0.220000 

7.816345 

28.137617 

6.562022 

8, 

0.240000 

8.398710 

30.123139 

7.256241 

8, 

0.260000 

9.022239 

32.253948 

7.984031 

9, 

0.280000 

9.689845 

34.531365 

8.764623 

10, 

0.300000 

10.404493 

36.958783 

9.612431 

11, 

0.320000 

11.169226 

39.540619 

10.539757 

12, 

0.340000 

11.987179 

42.281684 

11.558226 

13, 

0.360000 

12.861587 

45.186784 

12.679589 

14, 

0.380000 

13.795774 

48.260409 

13.916203 

16, 

0.400000 

14.793153 

51.506489 

15.281359 

17, 

0.420000 

15.857205 

54.928160 

16.789502 

19, 

0.440000 

16.991464 

58.527535 

18.456426 

20, 

0.460000 

18.199497 

62.305471 

20.299425 

22. 

0.480000 

19.484870 

66.261315 

22.337436 

24. 

0.500000 

20.851120 

70.392640 

24.591167 

27. 

0.520000 

22.301716 

74.694967 

27.083209 

29. 

0.540000 

23.840014 

79.161473 

29.838138 

32. 

0.560000 

25.469207 

83.782683 

32.882595 

35, 

0.580000 

27.192270 

88.546165 

36.245342 

38. 

0.600000 

29.011898 

93.436217 

39.957291 

42. 

0.620000 

30.930435 

98.433565 

44.051486 

46. 

0.640000 

32.949803 

103.515075 

48.563047 

51. 

0.660000 

35.071418 

108.653501 

53.529051 

56. 

0.680000 

37.296112 

113.817260 

58.988347 

61. 

0.700000 

39.624038 

118.970275 

64.981289 

67. 

0.720000 

42.054582 

124.071883 

71.549383 

74. 

0.740000 

44.586270 

129.076826 

78.734836 

81. 

0.760000 

47.216680 

133.935357 

86.579985 

89. 

0.780000 

49.942349 

138.593458 

95.126622 

97. 

0.800000 

52.758696 

142.993204 

104.415189 

107. 

0.820000 

55.659947 

147.073271 

114.483849 

117. 

0.840000 

58.639070 

150.769614 

125.367447 

127. 

0.860000 

61.687734 

154.016311 

137.096350 

139. 

0.880000 

64.796279 

156.746580 

149.695201 

152. 

0.900000 

67.953711 

158.893961 

163.181596 

165. 

0.920000 

71.147719 

160.393647 

177.564726 

179. 

0.940000 

74.364726 

161.183948 

192.844014 

195. 

0.960000 

77.589966 

161.207851 

209.007797 

211. 

0.980000 

80.807591 

160.414636 

226.032105 

228. 

1.000000 

84.000817 

158.761500 

243.879597 

245. 

The  eigenvalues  are: 


1 , real 


1 , imag 


1 

2 

3 

4 


0.297688 

•3.383566 

•5.224641 

•6.271663 


14.940623 

36.623520 

53.193862 

68.921562 


T' 

800633 

696872 

221464 

146670 

207361 

559011 

123600 

223709 

080123 

822837 

525750 

231595 

966489 

747909 

588976 

500790 

493725 

578184 

765064 

066074 

493960 

062693 

787630 

685656 

775326 

076985 

612886 

407291 

486553 

879180 

615857 

729442 

254902 

229203 

691119 

680972 

240279 

411292 

236438 

757634 

015482 

048340 

891278 

574915 

124169 

556935 

882718 

101267 

201259 

159066 

937695 


Unstable:  at  least  one  eigenvalue  is  greater  than  zero 
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